304 *
#include <deal.II/base/conditional_ostream.h>
305 *
#include <deal.II/base/function.h>
306 *
#include <deal.II/base/function_lib.h>
307 *
#include <deal.II/base/index_set.h>
308 *
#include <deal.II/base/parameter_acceptor.h>
309 *
#include <deal.II/base/quadrature_lib.h>
310 *
#include <deal.II/base/timer.h>
311 *
#include <deal.II/base/utilities.h>
313 *
#include <deal.II/distributed/grid_refinement.h>
314 *
#include <deal.II/distributed/tria.h>
316 *
#include <deal.II/dofs/dof_handler.h>
317 *
#include <deal.II/dofs/dof_renumbering.h>
318 *
#include <deal.II/dofs/dof_tools.h>
320 *
#include <deal.II/fe/fe_q.h>
321 *
#include <deal.II/fe/fe_system.h>
322 *
#include <deal.II/fe/fe_values.h>
324 *
#include <deal.II/grid/grid_generator.h>
325 *
#include <deal.II/grid/grid_tools.h>
327 *
#include <deal.II/lac/dynamic_sparsity_pattern.h>
328 *
#include <deal.II/lac/full_matrix.h>
329 *
#include <deal.II/lac/generic_linear_algebra.h>
330 *
#include <deal.II/lac/solver_richardson.h>
331 *
#include <deal.II/lac/vector.h>
333 *
#include <deal.II/numerics/data_out.h>
334 *
#include <deal.II/numerics/error_estimator.h>
335 *
#include <deal.II/numerics/vector_tools.h>
337 *
#include <
boost/archive/binary_iarchive.hpp>
338 *
#include <
boost/archive/binary_oarchive.hpp>
339 *
#include <
boost/iostreams/device/file.hpp>
340 *
#include <
boost/iostreams/filter/gzip.hpp>
341 *
#include <
boost/iostreams/filtering_stream.hpp>
344 *
#include <filesystem>
346 *
#include <iostream>
353 * This
namespace contains the implementation of the
parallel flow routing
354 * program. The program solves a model of water flow on terrain
using the
355 *
"downhill flow" method commonly used in hydrology. In
this approach, water
356 * flows from each grid
point in the steepest downhill direction determined by
357 * the terrain elevation, and the water flux distribution is computed by solving
358 * a steady-state flow conservation system in
parallel.
361 *
namespace ParallelFlowRouting
367 * The LA
namespace encapsulates the linear algebra library configuration.
368 * We use either
PETSc or Trilinos
for distributed sparse matrices and
369 * vectors, depending on what deal.II was compiled with.
PETSc is preferred
370 *
if available (and not
using complex
numbers), otherwise we fall back to
371 * Trilinos. This choice allows
for efficient
parallel solving of the large
372 * sparse systems that arise from discretizing the flow conservation
378 *
#
if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
379 *
!(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
381 *
# define USE_PETSC_LA
382 *
#elif defined(DEAL_II_WITH_TRILINOS)
385 *
# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
391 * We use block vectors and block sparse matrices to represent the distributed
392 * linear systems. Block structures allow us to organize
data logically and
393 * can improve performance in certain scenarios.
397 *
using MatrixType = LA::MPI::BlockSparseMatrix;
403 * <a name=
"parallel_flow_routing.cc-ColoradoTopography"></a>
404 * <h3>ColoradoTopography</h3>
408 * This
class represents the topographic elevation
data for Colorado, defined
409 * on the domain spanning 7 degrees longitude by 4 degrees latitude
410 * (from 109°W to 102°W and from 37°N to 41°N). The
class reads elevation
data
411 * from a gzip-compressed file, processes it to ensure it is suitable
for flow
412 * routing (by filling depressions), and then caches the processed
data for
413 * reuse in subsequent runs.
417 * The elevation
data is provided at a baseline resolution of 1800 meters per
418 * pixel. The n_refinements parameter allows scaling to finer resolutions by
419 * subdividing each baseline grid cell into 2^n_refinements smaller cells.
420 * This is useful
for studying the flow routing algorithm at different
425 * The
class inherits from the deal.II
Function class, so it can be used
427 * values on a finite element mesh. The function is evaluated on 3D mesh
428 * points by
first converting Cartesian coordinates (used on the mesh) to
429 * geographic coordinates (longitude and latitude in degrees), and then
430 * interpolating the stored elevation
data.
434 * A key feature of this class is the removal of local depressions (sinks)
435 * from the elevation
data. These depressions are problematic
for flow routing
436 * because water trapped in them does not flow anywhere. The class uses a
437 * priority-flood depression-filling algorithm to ensure that the resulting
438 * digital elevation model (DEM) has no interior depressions. This is done
439 * during construction in
parallel, with process 0 performing the computation
440 * and broadcasting the result to all other processes.
443 *
class ColoradoTopography : public
Function<3>
448 * Constructor that loads and processes the Colorado topography
data.
449 * - `mpi_communicator`: The
MPI communicator used
for parallel
450 * communication. Process 0 reads and processes the
data, then
451 * broadcasts it to all others.
452 * - `n_refinements`: The refinement
level. The mesh will have
453 * 7*2^n_refinements by 4*2^n_refinements cells.
456 *
ColoradoTopography(
const MPI_Comm mpi_communicator,
457 *
const unsigned int n_refinements);
461 * Return the elevation in meters at a given 3D
point on the Earth
's
462 * surface. The point is first converted from Cartesian coordinates (as used
463 * on the mesh) to geographic coordinates (longitude and latitude in
464 * degrees), and then the stored elevation data is interpolated.
468 * value(const Point<3> &p,
469 * const unsigned int /*component*/ = 0) const override;
474 * The InterpolatedUniformGridData object stores the actual elevation values
475 * on a uniform grid in longitude-latitude space.
478 * std::unique_ptr<const Functions::InterpolatedUniformGridData<2>> data;
482 * Count the number of local depressions in an elevation table. A depression
483 * is a grid point that is lower than all of its 8 neighbors (including
487 * static unsigned int
488 * count_depressions(const Table<2, double> &elevation_data);
492 * Fill (eliminate) local depressions in an elevation table using a
493 * priority-flood algorithm. After this function completes, no interior grid
494 * point will be lower than all of its 8 neighbors.
498 * fill_depressions(Table<2, double> &elevation_data);
505 * <a name="parallel_flow_routing.cc-ColoradoTopographyColoradoTopography"></a>
506 * <h3>ColoradoTopography::ColoradoTopography()</h3>
510 * The constructor is responsible for loading the Colorado topography data.
511 * The procedure it follows is:
515 * 1. On process 0 (the root MPI process), check if a cache file exists for
516 * the requested resolution. If not, read the original data from a
517 * gzip-compressed file, interpolate it to the desired resolution, fill
518 * depressions, and cache the result in a binary-serialized file. If the
519 * cache exists, read directly from it (this is much faster).
523 * 2. Broadcast the elevation data to all MPI processes.
527 * 3. Create an InterpolatedUniformGridData object that can efficiently
528 * evaluate elevation at arbitrary points via bilinear interpolation.
532 * The input data file is expected to be in ESRI ASCII raster format,
533 * compressed with gzip. The format includes header lines specifying the
534 * number of columns and rows, the corner coordinates, and the cell size.
537 * ColoradoTopography::ColoradoTopography(const MPI_Comm mpi_communicator,
538 * const unsigned int n_refinements)
540 * unsigned int n_latitudes = numbers::invalid_unsigned_int;
541 * unsigned int n_longitudes = numbers::invalid_unsigned_int;
542 * Point<2> lower_left_corner = numbers::signaling_nan<Point<2>>();
543 * double pixel_size = numbers::signaling_nan<double>();
545 * Table<2, double> elevation_data;
547 * const unsigned int root_process = 0;
548 * if (Utilities::MPI::this_mpi_process(mpi_communicator) == root_process)
550 * const std::string cache_filename =
551 * "colorado-topography-1800m.cache" + std::to_string(n_refinements);
553 * if (!std::filesystem::exists(cache_filename))
555 * const std::string original_data_filename =
556 * "colorado-topography-1800m.txt.gz";
557 * std::cout << " Reading original elevation data from file "
558 * << original_data_filename << std::endl;
560 * unsigned int n_original_latitudes;
561 * unsigned int n_original_longitudes;
562 * Point<2> original_lower_left_corner;
563 * double original_pixel_size;
565 * Table<2, double> original_elevation_data;
567 * boost::iostreams::filtering_istream in;
568 * in.push(boost::iostreams::basic_gzip_decompressor<>());
569 * in.push(boost::iostreams::file_source(original_data_filename));
574 * AssertThrow(word == "ncols",
576 * "The first line of the input file needs to start "
577 * "with the word 'ncols
', but starts with '" +
579 * in >> n_original_longitudes;
582 * AssertThrow(word == "nrows",
584 * "The second line of the input file needs to start "
585 * "with the word 'nrows
', but starts with '" +
587 * in >> n_original_latitudes;
590 * AssertThrow(word == "xllcorner",
592 * "The third line of the input file needs to start "
593 * "with the word 'xllcorner
', but starts with '" +
595 * in >> original_lower_left_corner[0];
598 * AssertThrow(word == "yllcorner",
600 * "The fourth line of the input file needs to start "
601 * "with the word 'yllcorner
', but starts with '" +
603 * in >> original_lower_left_corner[1];
607 * AssertThrow(word == "cellsize",
609 * "The fourth line of the input file needs to start "
610 * "with the word 'cellsize
', but starts with '" +
612 * in >> original_pixel_size;
614 * original_elevation_data.reinit(n_original_longitudes,
615 * n_original_latitudes);
619 * Data is provided in the input file as horizontal strips, with
620 * longitude marching fastest from west to east. The second
621 * coordinate is latitude, but the file marches north to south,
622 * which is the opposite of how we want it (we want a right-handed
623 * coordinate system), so we have to revert the order in which we
624 * insert things into the data table.
627 * for (unsigned int latitude_index = 0;
628 * latitude_index < n_original_latitudes;
630 * for (unsigned int longitude_index = 0;
631 * longitude_index < n_original_longitudes;
639 * original_elevation_data(longitude_index,
640 * n_original_latitudes -
641 * latitude_index - 1) = elevation;
647 * "Could not read all expected data points "
648 * "from the file <" +
649 * original_data_filename + ">!"));
653 * const Functions::InterpolatedUniformGridData<2>
654 * original_elevation_field(
655 * std::array<std::pair<double, double>, 2>{
656 * {std::make_pair(original_lower_left_corner[0],
657 * original_lower_left_corner[0] +
658 * (n_original_longitudes - 1) *
659 * original_pixel_size),
660 * std::make_pair(original_lower_left_corner[1],
661 * original_lower_left_corner[1] +
662 * (n_original_latitudes - 1) *
663 * original_pixel_size)}},
664 * std::array<unsigned int, 2>{
665 * {n_original_longitudes - 1, n_original_latitudes - 1}},
666 * std::move(original_elevation_data));
670 * The model that we just read is provided on a domain that is
671 * slightly larger than what we actually need, and on a mesh that
672 * does not align with the vertices we will create later on. It may
673 * also contain local depressions that prevent us from performing
674 * useful flow routing, and even if it doesn't, the interpolation
675 * onto a concrete mesh that has different vertices will create a
676 * model with local depressions.
680 * To avoid
this, we take the following steps where
'r' is a
681 * parameter that controls the resolution of the mesh we will
683 * * We
interpolate things onto a mesh that has 7*2^r x 4*2^r cells
684 * (i.e., 7*2^r+1 x 4*2^r+1 points) and that has the exact right
685 * extents. This makes sense because below we will
686 * start all computations on a 7x4 mesh, given that Colorado spans
687 * 7 by 4 degrees on the surface of the Earth.
688 * * We then perform depression filling by lifting points that are
689 * lower than all of their neighbors above the lowest of its
690 * neighbors. The result is a digital elevation model without
694 *
const unsigned int n_subdivisions = (1 << n_refinements);
695 *
n_longitudes = 7 * n_subdivisions + 1;
696 *
n_latitudes = 4 * n_subdivisions + 1;
697 *
pixel_size = 1. / n_subdivisions;
698 *
lower_left_corner = {-109., 37.};
700 *
elevation_data.reinit(n_longitudes, n_latitudes);
701 *
for (
unsigned int latitude_index = 0; latitude_index < n_latitudes;
703 *
for (
unsigned int longitude_index = 0;
704 *
longitude_index < n_longitudes;
707 *
const double longitude =
708 *
lower_left_corner[0] + longitude_index * pixel_size;
709 *
const double latitude =
710 *
lower_left_corner[1] + latitude_index * pixel_size;
711 *
elevation_data(longitude_index, latitude_index) =
712 *
original_elevation_field.value(
718 * Now we need to fix up depressions in the *interior* of the table.
721 *
std::cout <<
" Filling in the "
722 *
<< count_depressions(elevation_data)
723 *
<<
" depressions in the elevation model" << std::endl;
724 *
fill_depressions(elevation_data);
728 * Check that we have no depressions left:
731 *
Assert(count_depressions(elevation_data) == 0, ExcInternalError());
735 * Write the elevation
data to the cache file
using binary
739 *
std::cout <<
" Writing " << elevation_data.size()[0] <<
" x "
740 *
<< elevation_data.size()[1]
741 *
<<
" elevation points to cache file " << cache_filename
743 *
boost::iostreams::filtering_ostream out;
744 *
out.push(boost::iostreams::basic_gzip_compressor<>());
745 *
out.push(boost::iostreams::file_sink(cache_filename));
746 *
boost::archive::binary_oarchive oa(out);
747 *
oa << n_longitudes << n_latitudes << lower_left_corner << pixel_size
754 * If we did find a cache file read by a previous
run of the
756 * read the elevation
data using binary deserialization
759 *
std::cout <<
" Reading elevation data from cache file "
760 *
<< cache_filename << std::endl;
762 *
boost::iostreams::filtering_istream in;
763 *
in.push(boost::iostreams::basic_gzip_decompressor<>());
764 *
in.push(boost::iostreams::file_source(cache_filename));
765 *
boost::archive::binary_iarchive ia(in);
766 *
ia >> n_longitudes >> n_latitudes >> lower_left_corner >>
767 *
pixel_size >> elevation_data;
769 *
std::cout <<
" Read " << elevation_data.size()[0] <<
" x "
770 *
<< elevation_data.size()[1]
771 *
<<
" elevation points from cache file" << std::endl;
777 * Finally, distribute the
data created on process 0 to everyone
else and
778 * create a function
object that can be used to evaluate the elevation at
792 *
elevation_data.replicate_across_communicator(mpi_communicator,
795 *
data = std::make_unique<const Functions::InterpolatedUniformGridData<2>>(
796 *
std::array<std::pair<double, double>, 2>{
797 *
{std::make_pair(lower_left_corner[0],
798 *
lower_left_corner[0] + (n_longitudes - 1) * pixel_size),
799 *
std::make_pair(lower_left_corner[1],
800 *
lower_left_corner[1] +
801 *
(n_latitudes - 1) * pixel_size)}},
802 *
std::array<unsigned int, 2>{{n_longitudes - 1, n_latitudes - 1}},
803 *
std::move(elevation_data));
810 * <a name=
"parallel_flow_routing.cc-ColoradoTopographyvalue"></a>
811 * <h3>ColoradoTopography::value()</h3>
815 * This function evaluates the elevation at a given 3D
point. Since the mesh
816 * is embedded in 3D on the surface of the Earth (as a sphere of
radius 6371
817 * km), the input
point p is given in Cartesian coordinates. We convert these
818 * to geographic coordinates (longitude and latitude in degrees) and then look
819 * up the elevation in the stored
data table
using bilinear interpolation.
823 * The conversion from Cartesian to geographic coordinates uses standard
825 * - Longitude (x-y plane angle):
atan2(y, x) * 360 / (2π)
826 * - Latitude (z-radial
angle):
atan2(z,
sqrt(x² + y²)) * 360 / (2π)
835 * First pull back p to longitude/latitude, expressed in degrees
841 *
std::
sqrt(p[0] * p[0] + p[1] * p[1])) *
844 *
return data->value(p_long_lat);
851 * <a name=
"parallel_flow_routing.cc-ColoradoTopographycount_depressions"></a>
852 * <h3>ColoradoTopography::count_depressions()</h3>
856 * This
static helper function counts the number of local depressions in an
857 * elevation table.
A local depression is a grid
point in the interior of the
858 * domain that is lower than all of its 8 neighbors (i.e., all
859 * immediate neighbors including diagonals).
863 * We only
check interior points (excluding the boundary of the domain)
864 * because boundary points can have lower neighbors (water can flow out of the
865 * domain at the boundary at these points).
869 *
ColoradoTopography::count_depressions(
const Table<2, double> &elevation_data)
871 *
const unsigned int n_longitudes = elevation_data.size()[0];
872 *
const unsigned int n_latitudes = elevation_data.size()[1];
874 *
unsigned int n_depressions = 0;
875 *
for (
unsigned int x = 1; x < n_longitudes - 1; ++x)
876 *
for (
unsigned int y = 1; y < n_latitudes - 1; ++y)
878 *
const double elevation = elevation_data(x, y);
879 *
double min_neighbor_elevation = std::numeric_limits<double>::max();
880 *
for (
int i = -1; i <= +1; ++i)
881 *
for (
int j = -1; j <= +1; ++j)
882 *
if (!(i == 0 && j == 0))
883 *
min_neighbor_elevation =
std::min(min_neighbor_elevation,
884 *
elevation_data(x + i, y + j));
885 *
if (min_neighbor_elevation >= elevation)
888 *
return n_depressions;
895 * <a name=
"parallel_flow_routing.cc-ColoradoTopographyfill_depressions"></a>
896 * <h3>ColoradoTopography::fill_depressions()</h3>
900 * This function removes local depressions (sinks) from a gridded elevation
901 * model
using a priority-flood algorithm. In topographic
data, depressions
902 * are local minima that are lower than all their neighbors; they are
903 * problematic
for flow routing because they trap water that would
904 * otherwise flow downhill.
908 * The key idea is to work
"inward" from the boundaries: we start by marking
909 * all border cells as processed and placing them in a priority queue
sorted
910 * by elevation (lowest
first). Then, we repeatedly
extract the lowest cell
911 * from the queue and examine its unprocessed 8-connected neighbors
912 * (including diagonals). For each unprocessed neighbor, we compute its
913 * filled elevation as the maximum of its current elevation and the parent
914 * cell
's elevation. To avoid creating perfectly flat plateaus that would
915 * be ambiguous for flow routing, we add a small deterministic increment
916 * (derived from the cell's grid indices). We then mark the neighbor as
917 * processed, update its elevation in-place, and add it to the queue. This
918 * continues until the queue is empty.
922 * The result is that no interior cell is lower than all of its neighbors.
923 * The algorithm runs in
O(N log N) time where
N is the number of grid cells,
924 * and the deterministic increment ensures reproducibility across runs.
928 * @note In hindsight, the choice of a
random increment between
zero and
929 * 0.1 may have been a bit large -- we really just want to avoid flat
930 * plateaus, so an increment on the order of 0.01 or even smaller would have
931 * been sufficient. The current choice may create some small artificial
932 * slopes that might add up to too much elevation change across the domain.
933 * But,
this is the
value used
for the experiments in the accompanying paper,
934 * so we keep it as is
for now.
947 *
return elev > other.elev;
951 *
const unsigned int n_rows = elevation_data.size()[0];
952 *
const unsigned int n_cols = elevation_data.size()[1];
954 *
processed.fill(
false);
958 *
A priority queue that
always gives us the lowest elevation node that
959 * has not yet been processed. std::priority_queue is a rarely used
960 *
data structure that is a heap-based implementation of a priority queue.
961 * You can find its description at
965 *
std::priority_queue<Node,
967 *
std::greater<Node>>
968 *
currently_active_nodes;
974 * Push all border nodes into the priority queue
977 *
for (
unsigned int i = 0; i < n_rows; ++i)
979 *
currently_active_nodes.push(Node{i, 0, elevation_data[i][0]});
980 *
currently_active_nodes.push(
981 *
Node{i, n_cols - 1, elevation_data[i][n_cols - 1]});
982 *
processed[i][0] = processed[i][n_cols - 1] =
true;
984 *
for (
unsigned int j = 0; j < n_cols; ++j)
986 *
currently_active_nodes.push(Node{0, j, elevation_data[0][j]});
987 *
currently_active_nodes.push(
988 *
Node{n_rows - 1, j, elevation_data[n_rows - 1][j]});
989 *
processed[0][j] = processed[n_rows - 1][j] =
true;
994 * Directions
for 8 neighbors
997 *
const int dx[8] = {-1, -1, -1, 0, 0, 1, 1, 1};
998 *
const int dy[8] = {-1, 0, 1, -1, 1, -1, 0, 1};
1002 * While there are nodes
for which we can still look
for neighbors:
1005 *
while (!currently_active_nodes.empty())
1007 *
const Node current_node =
1008 *
currently_active_nodes
1010 *
currently_active_nodes.pop();
1014 * Loop over the neighbors of the lowest
point in the queue, excluding
1015 * points that are beyond the boundary of the domain and also skipping
1016 * over the ones that have already been processed:
1019 *
for (
unsigned int k = 0; k < 8; ++k)
1020 *
if ((
static_cast<signed int>(current_node.x) + dx[k] >= 0) &&
1021 *
(
static_cast<signed int>(current_node.x) + dx[k] <
1022 *
static_cast<signed int>(n_rows)) &&
1023 *
(
static_cast<signed int>(current_node.y) + dy[k] >= 0) &&
1024 *
(
static_cast<signed int>(current_node.y) + dy[k] <
1025 *
static_cast<signed int>(n_cols)))
1028 *
const unsigned int neighbor_x = current_node.x +
dx[k];
1029 *
const unsigned int neighbor_y = current_node.y + dy[k];
1031 *
if (processed[neighbor_x][neighbor_y])
1036 * If the neighbor exists and has not been processed:
1037 * Add a
random increment between 0 and 0.1 meters to ensure
1038 * that we
do not create perfectly flat plateaus. The increment is
1039 *
random but deterministic to ensure reproducibility across runs.
1042 *
const double new_elevation =
1043 *
std::max(elevation_data[neighbor_x][neighbor_y],
1044 *
current_node.elev +
1045 *
std::uniform_real_distribution<>(0, 0.1)(rng));
1046 *
elevation_data[neighbor_x][neighbor_y] = new_elevation;
1047 *
processed[neighbor_x][neighbor_y] =
true;
1051 * Push that neighbor to the queue:
1054 *
currently_active_nodes.push(
1055 *
{
static_cast<unsigned int>(neighbor_x),
1056 *
static_cast<unsigned int>(neighbor_y),
1067 * <a name=
"parallel_flow_routing.cc-RainFallRate"></a>
1068 * <h3>RainFallRate</h3>
1072 * This is a simple
class that describes the rain fall rate on the domain.
1073 * In reality, the rain fall rate varies depending on location and climate
1074 * conditions, but
for this program we use a
constant value everywhere on
1075 * the domain. The rain fall rate is an important boundary condition
for
1076 * the flow routing problem: it represents the water that enters the system
1077 * through precipitation, and eventually either flows out through the
1078 * boundary or accumulates in local depressions.
1082 * The
value of 375 mm per year (approximately 15 inches per year) is
1083 * roughly representative of the rainfall in Colorado.
1086 *
template <
int spacedim>
1087 *
class RainFallRate :
public Function<spacedim>
1092 *
const unsigned int component = 0)
const override;
1096 *
template <
int spacedim>
1099 *
const unsigned int )
const
1147 *
static constexpr int dim = 2;
1148 *
static constexpr int spacedim = 3;
1150 *
ParallelFlowRouter();
1163 *
interpolate_initial_elevation();
1166 *
sort_dofs_high_to_low();
1169 *
compute_local_flow_routing();
1172 *
assemble_system();
1175 *
assemble_matrix_free_operators();
1181 *
check_conservation_for_waterflow_system(
const VectorType &solution);
1188 *
unsigned int n_refinements;
1189 *
bool generate_graphical_output;
1196 *
std::vector<IndexSet> locally_owned_partitioning;
1197 *
std::vector<IndexSet> locally_relevant_partitioning;
1199 *
IndexSet locally_relevant_water_dofs;
1201 *
std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
1202 *
local_flow_routing;
1209 *
class FlowRoutingMatrix;
1210 *
std::unique_ptr<const FlowRoutingMatrix> flow_routing_matrix;
1212 *
class FlowRoutingPreconditioner;
1213 *
std::unique_ptr<const FlowRoutingPreconditioner>
1214 *
flow_routing_preconditioner;
1216 *
class IplusminusXMatrixBase;
1218 *
class IplusXMatrix;
1219 *
std::unique_ptr<const IplusXMatrix> I_plus_X_matrix;
1221 *
class IminusXMatrix;
1222 *
std::unique_ptr<const IminusXMatrix> I_minus_X_matrix;
1229 *
ParallelFlowRouter::ParallelFlowRouter()
1231 *
, mpi_communicator(MPI_COMM_WORLD)
1232 *
, n_refinements(9)
1234 *
, triangulation(mpi_communicator,
1238 *
, dof_handler(triangulation)
1239 *
, pcout(std::cout,
1241 *
, computing_timer(mpi_communicator,
1246 *
add_parameter(
"Number of refinements",
1248 *
"The number of global refinements to perform.");
1249 *
add_parameter(
"Generate graphical output",
1250 *
generate_graphical_output,
1251 *
"Whether to generate graphical output files.");
1258 * <a name=
"parallel_flow_routing.cc-ParallelFlowRoutermake_grid"></a>
1259 * <h3>ParallelFlowRouter::make_grid()</h3>
1263 * This function creates the mesh that discretizes the Colorado topography
1264 * domain. Rather than
using a flat Cartesian coordinate system, the mesh is
1265 * mapped onto the surface of the Earth (a sphere of
radius 6371 km), so that
1266 * distances and areas computed on the mesh are realistic.
1270 * The function starts by creating a rectangular mesh that spans Colorado
's
1271 * geographic extent (7 degrees longitude × 4 degrees latitude). After
1272 * global refinements, this mesh is then transformed from longitude-latitude
1273 * coordinates to 3D Cartesian coordinates on the Earth's surface via a
1274 * cylindrical projection. This ensures that when we compute
gradients or
1275 * areas on the mesh, they respect the curvature of the Earth, which is
1276 * essential
for correct flow routing on real terrain.
1280 *
ParallelFlowRouter::make_grid()
1283 *
pcout <<
"Making grid... " << std::endl;
1291 *
triangulation.refine_global(n_refinements);
1295 *
const Point<2> p_long_lat(p_long_lat_degrees[0] / 360 *
1297 *
p_long_lat_degrees[1] / 360 *
1299 *
const double R = 6371000;
1308 *
pcout <<
" Number of cells: " << triangulation.n_global_active_cells()
1312 *
for (
const auto &cell : triangulation.active_cell_iterators())
1313 *
area += cell->measure();
1314 *
pcout <<
" Area of the domain: " << area <<
"m^2" << std::endl;
1321 * <a name=
"parallel_flow_routing.cc-ParallelFlowRoutersetup_system"></a>
1322 * <h3>ParallelFlowRouter::setup_system()</h3>
1326 * The next function also is not something that is
new in any particular
1327 * way. Conceptually, all we have to
do is
set up block vectors and matrices
1328 *
for the linear systems we want to solve. This really is quite
1329 * straightforward with the only complication that we have to account
for
1330 * the fact that we are working in a
parallel program where we have to
1331 * keep track which process owns which degrees of freedom.
1335 * This function does the basic
set-up. In following
functions, we will
1336 * re-enumerate the degrees of freedom in a way that makes the flow routing
1337 *
matrix have a nice triangular structure, and we will also
set up the
1338 *
matrix-free operators that we will use to compare the
matrix-based and
1343 *
ParallelFlowRouter::setup_dofs()
1346 *
pcout <<
"Setting up system... " << std::endl;
1348 *
dof_handler.distribute_dofs(fe);
1351 *
const std::vector<types::global_dof_index> dofs_per_block =
1358 *
const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
1359 *
locally_owned_partitioning = {
1360 *
locally_owned_dofs.
get_view(0, n_elevation_dofs),
1361 *
locally_owned_dofs.
get_view(n_elevation_dofs,
1362 *
n_elevation_dofs + n_waterflow_rate_dofs)};
1366 *
locally_relevant_dofs =
1368 *
locally_relevant_partitioning = {
1369 *
locally_relevant_dofs.get_view(0, n_elevation_dofs),
1370 *
locally_relevant_dofs.get_view(n_elevation_dofs,
1371 *
n_elevation_dofs + n_waterflow_rate_dofs)};
1373 *
IndexSet all_elevation_dofs(dof_handler.n_dofs());
1374 *
all_elevation_dofs.add_range(0, n_elevation_dofs);
1375 *
locally_relevant_water_dofs = locally_relevant_dofs;
1376 *
locally_relevant_water_dofs.subtract_set(all_elevation_dofs);
1378 *
locally_relevant_solution.reinit(locally_owned_partitioning,
1379 *
locally_relevant_partitioning,
1380 *
mpi_communicator);
1381 *
locally_relevant_solution_dot.reinit(locally_owned_partitioning,
1382 *
locally_relevant_partitioning,
1383 *
mpi_communicator);
1384 *
system_rhs.reinit(locally_owned_partitioning, mpi_communicator);
1386 *
pcout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
1387 *
<<
" (elevation: " << n_elevation_dofs
1388 *
<<
", waterflow: " << n_waterflow_rate_dofs <<
')' << std::endl;
1389 *
const std::vector<types::global_dof_index> water_dofs_per_process =
1391 *
locally_owned_partitioning[1].n_elements());
1392 *
pcout <<
" Number of waterflow degrees of freedom per process: "
1393 *
<< std::accumulate(water_dofs_per_process.begin(),
1394 *
water_dofs_per_process.end(),
1397 *
<<
" (average) x "
1399 *
<<
" (number of processes)" << std::endl;
1406 * <a name=
"parallel_flow_routing.cc-ParallelFlowRouterinterpolate_initial_elevation"></a>
1407 * <h3>ParallelFlowRouter::interpolate_initial_elevation()</h3>
1411 * The following function then interpolates the
initial elevation onto the
1412 * mesh. We need
this as
initial conditions
for the elevation variable.
1416 * The way
this function works is that given a
scalar function object (derived
1417 * from the
Function class), we
first create a function
object that covers
1418 * all solution variables, returns the elevation in
one vector
component
1419 * (specifically, in vector
component zero) and zeros in all others. This
1420 *
"extension" of a scalar to a vector function is done by the
1422 * extended function
object to interpolate these initial conditions onto
1423 * all degrees of freedom, which correctly sets the initial elevation
1424 * variables to their initial values and sets the water flow rate to
1429 *
ParallelFlowRouter::interpolate_initial_elevation()
1432 *
"Initial conditions: interpolate elevation");
1433 *
pcout <<
"Interpolating elevation... " << std::endl;
1435 *
const ColoradoTopography colorado_topography(mpi_communicator,
1438 *
[&](
const Point<spacedim> &p) {
return colorado_topography.value(p); },
1442 *
VectorType interpolated_initial_condition(locally_owned_partitioning,
1446 *
interpolated_initial_condition);
1450 * The vector we have just interpolated into is a
"fully distributed
1451 * vector", i.e., every element is uniquely owned by
one of the
MPI
1452 * processes and these are the only ones we store on the current process.
1453 * On the other hand, we will also have to access
values for nodes that
1454 * are owned by other processes (
for example on ghost cells), so we
copy
1455 * the vector into
one that also has these ghost entries:
1458 *
locally_relevant_solution = interpolated_initial_condition;
1466 * <a name=
"parallel_flow_routing.cc-ParallelFlowRoutersort_dofs_high_to_low"></a>
1467 * <h3>ParallelFlowRouter::sort_dofs_high_to_low()</h3>
1471 * The key insight
for efficiently solving the flow routing problem is that
1472 * water
always flows downhill. This means we can process water flow
1473 * calculations in order from the highest elevation to the lowest. This
1474 * ordering is crucial
for the performance of the solver: it means that when
1475 * we compute the water flow at a given node, all of the upstream nodes from
1476 * which it receives water have already been processed. This is the
1477 * characteristic of a triangular system, and it allows us to solve the
1478 * system very efficiently
using a simple substitution method.
1482 * This function implements
this idea by renumbering the degrees of freedom
1483 * that represent the water flow rate such that they are ordered from highest
1484 * elevation to lowest. This doesn
't change the mathematical problem we're
1485 * solving, but it transforms the
matrix of the linear system into a lower
1486 * triangular
matrix when water flows are processed in
this order, making
1487 * the solver much more efficient.
1491 * The algorithm works by
first collecting all water degrees of freedom and
1492 * their corresponding elevations, sorting them from highest to lowest, and
1493 * then renumbering the DoFs accordingly. After renumbering, we update the
1494 *
index sets used
for parallel communication to reflect the
new ordering.
1498 * As the paper notes,
this reordering step is not actually necessary
if
1499 * you implement the per-process high-to-low flow routing algorithm in a
1500 *
matrix-free way, but it is convenient
for the
matrix-based solver because
1501 * in that
case, the preconditioner is a triangular solve which is exactly
1502 * what the Gauss-Seidel (SOR) method will perform and so we can use what
1503 *
PETSc offers to us without having to implement a custom preconditioner.
1507 *
ParallelFlowRouter::sort_dofs_high_to_low()
1510 *
"Initial conditions: sort DoFs high to low");
1511 *
pcout <<
"Sorting DoFs high to low... " << std::endl;
1513 *
std::map<types::global_dof_index, double> water_dof_index_to_elevation_map;
1514 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1515 *
if (cell->is_locally_owned())
1516 *
for (unsigned
int v = 0; v < cell->reference_cell().n_vertices(); ++v)
1519 *
cell->vertex_dof_index(v, 1);
1521 *
if (dof_handler.locally_owned_dofs().is_element(vertex_water_dof))
1524 *
cell->vertex_dof_index(v, 0);
1525 *
const double vertex_elevation =
1526 *
locally_relevant_solution(vertex_elevation_dof);
1527 *
water_dof_index_to_elevation_map[vertex_water_dof] =
1532 *
std::vector<std::pair<types::global_dof_index, double>>
1533 *
water_dof_index_to_elevation_list(
1534 *
water_dof_index_to_elevation_map.size());
1535 *
std::copy(water_dof_index_to_elevation_map.begin(),
1536 *
water_dof_index_to_elevation_map.end(),
1537 *
water_dof_index_to_elevation_list.begin());
1538 *
std::sort(water_dof_index_to_elevation_list.begin(),
1539 *
water_dof_index_to_elevation_list.end(),
1540 *
[](
const std::pair<types::global_dof_index, double> &a,
1541 *
const std::pair<types::global_dof_index, double> &b) {
1542 * return a.second > b.second;
1548 * low based on their elevation. We want to renumber the existing water
1549 * DoF indices so that they follow
this ordering. We will
do so
1550 * by making use of the fact that locally owned DoFs are numbered in
1551 * a
contiguous block, so we can start numbering all locally owned
1555 *
Assert(locally_owned_partitioning[0].is_contiguous(), ExcInternalError());
1556 *
Assert(locally_owned_partitioning[1].is_contiguous(), ExcInternalError());
1558 *
water_dof_index_to_elevation_list.size());
1559 *
std::map<types::global_dof_index, types::global_dof_index>
1560 *
old_to_new_water_indices;
1561 *
if (water_dof_index_to_elevation_list.size() >
1565 *
locally_owned_partitioning[0].size() +
1566 *
*locally_owned_partitioning[1].begin();
1567 *
for (
unsigned int i = 0; i < water_dof_index_to_elevation_list.size();
1569 *
old_to_new_water_indices[water_dof_index_to_elevation_list[i].
first] =
1570 *
first_water_dof_index + i;
1575 * Now we can
do the renumbering:
1578 *
std::vector<types::global_dof_index> new_dof_numbers;
1579 *
new_dof_numbers.reserve(dof_handler.n_locally_owned_dofs());
1583 * Do not re-enumerate the elevation DoFs at all.
1587 *
locally_owned_partitioning[0])
1588 *
new_dof_numbers.
push_back(elevation_dof);
1591 *
dof_handler.locally_owned_dofs())
1592 *
if (dof_index >= dof_handler.n_dofs() / 2)
1594 *
Assert(old_to_new_water_indices.find(dof_index) !=
1595 *
old_to_new_water_indices.
end(),
1596 *
ExcInternalError());
1597 *
new_dof_numbers.push_back(old_to_new_water_indices[dof_index]);
1599 *
AssertDimension(new_dof_numbers.size(), dof_handler.n_locally_owned_dofs());
1601 *
dof_handler.renumber_dofs(new_dof_numbers);
1605 * Rebuild
index sets after renumbering
1608 *
locally_relevant_dofs =
1610 *
const std::vector<types::global_dof_index> dofs_per_block =
1613 *
IndexSet all_elevation_dofs(dof_handler.n_dofs());
1614 *
all_elevation_dofs.add_range(0, n_elevation_dofs);
1615 *
locally_relevant_water_dofs = locally_relevant_dofs;
1616 *
locally_relevant_water_dofs.subtract_set(all_elevation_dofs);
1617 *
locally_relevant_partitioning = {
1618 *
locally_relevant_dofs.get_view(0, n_elevation_dofs),
1619 *
locally_relevant_dofs.get_view(n_elevation_dofs, dof_handler.n_dofs())};
1622 *
<<
" Elevations range between "
1624 *
water_dof_index_to_elevation_list.back().second :
1625 *
std::numeric_limits<double>::
max(),
1629 *
water_dof_index_to_elevation_list.
size() > 0 ?
1630 *
water_dof_index_to_elevation_list.front().
second :
1631 *
-
std::numeric_limits<double>::
max(),
1633 *
<<
"m." <<
std::endl;
1640 * <a name=
"parallel_flow_routing.cc-ParallelFlowRoutercompute_local_flow_routing"></a>
1641 * <h3>ParallelFlowRouter::compute_local_flow_routing()</h3>
1645 * Next, we need to have a function that
for each degree of freedom (locally
1646 * owned or ghost) finds which downhill neighbor (
if any) it gives water to.
1647 * This corresponds to the D8 scheme used in many flow routing codes, where
1648 * each node gives water to exactly
one other node -- specifically, to the
1649 * neighbor in the direction of steepest descent. Because we consider the
1650 * four immediate neighbors of a node on a regular mesh plus the four
1651 *
diagonal neighbors (8 neighbors total, hence
"D8"), the neighbor with the
1652 * steepest downhill slope may not necessarily be the lowest-lying neighbor.
1656 * The algorithm works by examining all locally relevant cells and their
1657 * vertices. For each vertex, we compute the slope to all other vertices in
1658 * the same cell. If the slope is
negative (i.e., downhill), we
check if it
's
1659 * steeper than the previously found downhill direction, and if so, we record
1660 * it as the new destination for water from this vertex.
1664 * Since we're working in
parallel, we only have complete information about
1665 * the neighbors of locally owned vertices. For vertices on the boundaries
1667 * information with neighboring processes via ghost cells. This is done
using
1672 * The final result is stored in the local_flow_routing member variable,
1673 * which maps each locally relevant water degree of freedom to its downstream
1674 * neighbor (or to
numbers::invalid_dof_index if the vertex is at a boundary
1675 * or in a depression with no outlet).
1679 *
ParallelFlowRouter::compute_local_flow_routing()
1682 *
pcout <<
"Computing local routing... " << std::endl;
1686 * First, create a map from each DoF to a pair (index,slope) of neighbors.
1687 * These slopes must necessarily be
negative because water only flows
1692 * We initialize
this map by looping over all locally relevant water DoFs
1693 * and setting the
value of the map to (index=invalid, slope=0).
1697 *
std::pair<types::global_dof_index, double>>
1698 *
water_dofs_to_steepest_downhill_neighbor_and_slope;
1700 *
water_dofs_to_steepest_downhill_neighbor_and_slope[i] = {
1705 * Then
loop over all locally owned and ghost cells, get the locations
1706 * and elevations of the four vertices of each cell, along with the
1707 * indices of the water DoF.
1710 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1711 *
if (cell->is_locally_owned() || cell->is_ghost())
1715 *
std::array<types::global_dof_index, 4> vertex_water_dof_indices;
1716 *
std::array<Point<spacedim>, 4> vertex_locations;
1717 *
std::array<double, 4> vertex_elevations;
1721 *
vertex_water_dof_indices[v] = cell->vertex_dof_index(v, 1);
1722 *
vertex_locations[v] = cell->vertex(v);
1725 *
cell->vertex_dof_index(v, 0);
1726 *
vertex_elevations[v] =
1727 *
locally_relevant_solution(vertex_elevation_dof_index);
1732 * Next, determine the slope from each vertex to each of the other
1734 * the previously most downhill slope we have encountered
for a DoF,
1735 * then use
this as the direction in which
this DoF will give water.
1736 * (Of course, it is possible that we encounter a steeper downhill
1737 * direction next on
this cell, or on another cell; in that
case,
1738 * we will simply overwrite what we determine here.)
1744 *
if (vertex_elevations[
w] < vertex_elevations[v])
1748 * Compute the slope between two vertices. The slope is the
1749 * elevation difference divided by the distance. Since the
1750 * mesh is mapped onto the surface of the Earth (a sphere),
1752 * distance between two points would be along the surface of
1753 * the Earth. This would correctly account
for the Earth
's
1754 * curvature and give us the true slope on the terrain. But
1756 * also difficult to compute, so we simply use the
1757 * straight-line distance between points (through the
1758 * Earth), which is a good approximation
for small
1762 *
const double slope =
1763 *
(vertex_elevations[
w] - vertex_elevations[v]) /
1764 *
((vertex_locations[v] - vertex_locations[
w]).
norm());
1765 *
Assert(slope < 0, ExcInternalError());
1767 *
water_dofs_to_steepest_downhill_neighbor_and_slope
1768 *
[vertex_water_dof_indices[v]]
1770 *
water_dofs_to_steepest_downhill_neighbor_and_slope
1771 *
[vertex_water_dof_indices[v]] = {
1772 *
vertex_water_dof_indices[
w], slope};
1778 * At
this point, we no longer care about slopes because we have considered
1779 * all neighbors of all nodes and no longer need to compare slopes between
1780 * nodes and neighbors. So
reduce the map to a smaller
one that only
1781 * contains
for each DoF who it gives water to.
1785 * Secondly, we have worked on all locally relevant DoFs up to
this point.
1786 * For all locally active DoFs (locally owned plus the ones on the interface
1787 * to ghost cells), we have considered all neighboring cells and so we can
1788 * be certain that we have their
downstream neighbors right. But
for the
1789 * nodes on the far side of the ghost cells (adjacent to artificial cells),
1790 * we have not seen all neighbor nodes, and so might have gotten wrong who
1791 * they give water to. As a consequence, we exclude those DoFs from the
1792 * reduced list that are not locally owned and instead obtain their
1793 * information via a ghost exchange. (We could exclude only the ones that
1794 * are not locally *active*, but that doesn
't buy us anything and the test
1795 * for locally owned is cheaper because that's a
contiguous set.)
1798 *
std::map<types::global_dof_index, types::global_dof_index>
1799 *
water_dofs_to_steepest_downhill_neighbor;
1800 *
for (
const auto &[source_index, dest_index_and_slope] :
1801 *
water_dofs_to_steepest_downhill_neighbor_and_slope)
1802 *
if (dof_handler.locally_owned_dofs().is_element(source_index))
1803 *
water_dofs_to_steepest_downhill_neighbor.
insert(
1804 *
{source_index, dest_index_and_slope.
first});
1807 *
dof_handler.locally_owned_dofs().n_elements() / 2);
1809 *
using CellLocalData =
1810 *
std::map<types::global_dof_index, types::global_dof_index>;
1814 * Pack up the locally owned water dof
index entries in the map
1815 * above
for the current cell:
1818 *
const auto pack_function =
1819 *
[
this, &water_dofs_to_steepest_downhill_neighbor](
1821 *
Assert(cell->is_locally_owned(), ExcInternalError());
1823 *
CellLocalData cell_local_data;
1827 *
cell->vertex_dof_index(v, 1);
1828 *
if (dof_handler.locally_owned_dofs().is_element(
1829 *
vertex_water_dof_index))
1831 *
Assert(water_dofs_to_steepest_downhill_neighbor.find(
1832 *
vertex_water_dof_index) !=
1833 *
water_dofs_to_steepest_downhill_neighbor.end(),
1834 *
ExcInternalError());
1835 *
cell_local_data.insert({vertex_water_dof_index,
1836 *
water_dofs_to_steepest_downhill_neighbor
1837 *
[vertex_water_dof_index]});
1840 *
return cell_local_data;
1845 * Unpack what the other processes have sent
for the current cell (which
1846 * is a ghost cell here). Because these were locally owned on the other
1847 * cell, they are necessarily not locally owned but locally relevant
1848 * here, and we assert that.
1852 * We will ultimately only care about flow from
one to another node
if
1853 * at least
one of them is locally active. We already know that the
1854 * source
index is not locally active, so we discard entries that have
1855 * a destination that is not locally active either. (You
'd think we
1856 * could have filtered this out in the pack_function above already, but
1857 * what we pack up on a cell may be sent to multiple processes that have
1858 * this cell as a ghost cell, and while a destination index may not be
1859 * locally relevant on one process, it may be on another.) We keep
1860 * the ones where the destination is not set, which indicates that
1861 * the DoF is at the boundary or in a depression without outlet.
1864 * const auto unpack_function =
1865 * [this, &water_dofs_to_steepest_downhill_neighbor](
1866 * const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
1867 * const CellLocalData &cell_local_data) {
1868 * Assert(cell->is_ghost(), ExcInternalError());
1870 * for (const auto &[source_index, dest_index] : cell_local_data)
1872 * Assert(dof_handler.locally_owned_dofs().is_element(source_index) ==
1874 * ExcInternalError());
1875 * Assert(locally_relevant_dofs.is_element(source_index) == true,
1876 * ExcInternalError());
1877 * if ((dest_index == numbers::invalid_dof_index) ||
1878 * locally_relevant_dofs.is_element(dest_index))
1879 * water_dofs_to_steepest_downhill_neighbor.insert(
1880 * {source_index, dest_index});
1884 * GridTools::exchange_cell_data_to_ghosts<CellLocalData>(dof_handler,
1890 * At this point, we should have gotten information about all locally
1891 * relevant water DoFs where they send their water (if anywhere),
1892 * excluding not locally owned ones that sent water to not locally
1893 * relevant ones -- these are at the outer fringes of the ghost layer
1894 * sending water further afield. In other words, we need to have
1895 * information about all locally active ones and at least some of
1896 * the locally relevant ones. We can check that this is the case:
1901 * Assert(water_dofs_to_steepest_downhill_neighbor.size() <=
1902 * locally_relevant_water_dofs.n_elements(),
1903 * ExcInternalError());
1904 * for (const auto &[src, dst] : water_dofs_to_steepest_downhill_neighbor)
1905 * Assert(locally_relevant_water_dofs.is_element(src),
1906 * ExcInternalError());
1908 * const types::global_dof_index n_elevation_dofs =
1909 * dof_handler.n_dofs() / 2;
1910 * IndexSet all_elevation_dofs(dof_handler.n_dofs());
1911 * all_elevation_dofs.add_range(0, n_elevation_dofs);
1912 * IndexSet locally_active_water_dofs =
1913 * DoFTools::extract_locally_active_dofs(dof_handler);
1914 * locally_active_water_dofs.subtract_set(all_elevation_dofs);
1915 * for (const auto &locally_active_index : locally_active_water_dofs)
1916 * Assert(water_dofs_to_steepest_downhill_neighbor.find(
1917 * locally_active_index) !=
1918 * water_dofs_to_steepest_downhill_neighbor.end(),
1919 * ExcInternalError());
1925 * Up to this point, it was useful to work with a std::map, but ultimately
1926 * we want a faster representation. So convert things into a std::vector
1930 * local_flow_routing = {water_dofs_to_steepest_downhill_neighbor.begin(),
1931 * water_dofs_to_steepest_downhill_neighbor.end()};
1935 * Finally, we can check that the only depressions we have should be
1936 * on the boundary of the domain. Recall that we marked depressions
1937 * (i.e., nodes that don't give water to any lower-lying neighbor)
1938 * in the src->dst relationships by invalid
'dst' values. This
1939 * should only be the
case for 'src' nodes that are on the
1940 * boundary, and we can check that:
1948 *
for (
const auto &[src, dst] : local_flow_routing)
1950 *
Assert(boundary_nodes.is_element(src),
1951 *
ExcMessage(
"Found an interior depression in the DEM."));
1960 * <a name=
"parallel_flow_routing.cc-ParallelFlowRouterassemble_system"></a>
1961 * <h3>ParallelFlowRouter::assemble_system()</h3>
1965 * This function assembles the linear system that describes the steady-state
1966 * water flow on the landscape. The system is based on the principle of mass
1967 * conservation: at each
point, the water flowing out must
equal the water
1968 * flowing in (from rain and from upstream neighbors) minus any water that
1973 * The discretized system has the form:
1974 * w_i = r_i + sum_{j: j flows to i} w_j
1975 * where w_i is the water flow rate at node i, r_i is the rainfall at node i,
1976 * and the
sum is over all upstream nodes j that flow into node i.
1980 * In
matrix form,
this becomes:
1982 * where F is the flow routing matrix (defined by the local_flow_routing
1983 *
data), I is the identity
matrix,
w is the vector of water flow rates, and r
1984 * is the rainfall vector.
1988 * Since we
've renumbered the DoFs so that water flows from higher to lower
1989 * elevations, the matrix (I - F) is lower triangular, making it easy to
1994 * Because we have chosen to work with a 2x2 block system where the first
1995 * block corresponds to elevation and the second block corresponds to water
1996 * flow rate, the sparsity patterns and matrix assembly are a bit more
1997 * complicated than in a standard finite element code, but the underlying
1998 * principles are the same. We just have to translate indices correctly and
1999 * make sure to fill the right blocks of the matrix.
2003 * ParallelFlowRouter::assemble_system()
2005 * TimerOutput::Scope t(computing_timer, "Solver 1: Assemble system");
2006 * pcout << "Assembling linear system... " << std::endl;
2008 * BlockDynamicSparsityPattern dsp(locally_relevant_partitioning);
2009 * for (const types::global_dof_index water_dof_within_block_1 :
2010 * locally_owned_partitioning[1])
2011 * dsp.block(1, 1).add(water_dof_within_block_1, water_dof_within_block_1);
2013 * for (const auto &[src, dst] : local_flow_routing)
2014 * if (dof_handler.locally_owned_dofs().is_element(dst))
2015 * dsp.add(dst, src);
2016 * SparsityTools::distribute_sparsity_pattern(dsp,
2017 * dof_handler.locally_owned_dofs(),
2019 * locally_relevant_dofs);
2023 * Now fill matrix accordingly
2026 * system_matrix.reinit(locally_owned_partitioning, dsp, mpi_communicator);
2027 * for (const types::global_dof_index water_dof_within_block_1 :
2028 * locally_owned_partitioning[1])
2029 * system_matrix.block(1, 1).set(water_dof_within_block_1,
2030 * water_dof_within_block_1,
2031 * 1.); // 1s on the diagonal
2032 * for (const auto &[water_dof, lowest_neighbor] : local_flow_routing)
2033 * if (dof_handler.locally_owned_dofs().is_element(lowest_neighbor))
2034 * system_matrix.set(lowest_neighbor,
2036 * -1); // -1s for flow routing
2037 * system_matrix.compress(VectorOperation::insert);
2041 * Then also fill rhs vector:
2044 * const RainFallRate<spacedim> rainfall_rate;
2045 * for (const auto &cell : dof_handler.active_cell_iterators())
2046 * if (cell->is_locally_owned())
2047 * for (unsigned int v = 0; v < cell->reference_cell().n_vertices(); ++v)
2049 * const types::global_dof_index vertex_water_dof =
2050 * cell->vertex_dof_index(v, 1);
2051 * system_rhs(vertex_water_dof) +=
2052 * rainfall_rate.value(cell->vertex(v)) * cell->measure() /
2053 * cell->n_vertices();
2055 * system_rhs.compress(VectorOperation::add);
2062 * <a name="parallel_flow_routing.cc-ParallelFlowRouterFlowRoutingMatrix"></a>
2063 * <h3>ParallelFlowRouter::FlowRoutingMatrix</h3>
2067 * This is a matrix-free operator class that represents the matrix A from the
2068 * discussion in assemble_system(). Rather than storing the matrix explicitly,
2069 * this class implements only the matrix-vector product (via the vmult()
2070 * function), computing the result on the fly from the local_flow_routing data
2075 * Recall that the matrix A encodes the flow routing: A has a -1 in position
2076 * (i, j) if water from node j flows to node i, and 0 elsewhere, plus a +1
2077 * on the diagonal. Each column j has at most one non-zero entry other
2078 * that the diagonal entry (since each node gives water to at most
2079 * one downhill neighbor). More precisely, for each src->dst pair in
2080 * local_flow_routing, we have a -1 in position (dst, src) of the matrix.
2084 * The vmult() function computes A*X = (I-F)*x = (I*x - F*x) by first
2085 * copying x to y (implementing I*x), then subtracting the contributions
2086 * from F. To compute F*x efficiently, we iterate over the src->dst pairs,
2087 * and for each one where 'dst
' is in the locally owned range, we add
2088 * -x[src] to y[dst].
2091 * class ParallelFlowRouter::FlowRoutingMatrix
2094 * FlowRoutingMatrix(
2095 * const IndexSet &locally_owned_water_dofs,
2096 * const IndexSet &locally_relevant_water_dofs,
2097 * const MPI_Comm mpi_communicator,
2098 * const unsigned int water_dofs_offset,
2099 * const std::vector<std::pair<types::global_dof_index,
2100 * types::global_dof_index>> &local_flow_routing)
2101 * : x_with_ghosts(locally_owned_water_dofs,
2102 * locally_relevant_water_dofs,
2104 * , my_local_flow_routing(local_flow_routing)
2108 * We got the map from DoFs to downhill neighbors in global DoF
2109 * indices, but we need them in indices relative to the second
2110 * vector block (or the (1,1) matrix block). So shift, unless
2111 * the destination DoF is invalid, indicating that this source
2112 * DoF has no outlet (because it's a depression in the DEM,
2113 * or because it
's at the boundary).
2116 * for (auto &[src, dst] : my_local_flow_routing)
2118 * src -= water_dofs_offset;
2119 * Assert(locally_relevant_water_dofs.is_element(src),
2120 * ExcInternalError());
2122 * if (dst != numbers::invalid_dof_index)
2124 * dst -= water_dofs_offset;
2125 * Assert(locally_relevant_water_dofs.is_element(dst),
2126 * ExcInternalError());
2132 * If one looks at how the vmult() function below is implemented,
2133 * one realizes that we only need those src->dst relationships
2134 * where 'dst
' is a valid DoF index and is in fact in the
2135 * locally owned range (it is the row index in the matrix,
2136 * and consequently that part of the output vector we fill
2137 * on the current process). To make this cheaper, we erase all others
2138 * at this point, so we don't have to
check any more there:
2141 *
Assert(locally_owned_water_dofs.is_contiguous(), ExcInternalError());
2143 *
std::remove_if(my_local_flow_routing.begin(),
2144 *
my_local_flow_routing.end(),
2145 *
[&locally_owned_water_dofs](
2148 * const types::global_dof_index dst = src_dst.second;
2149 * return ((dst == numbers::invalid_dof_index) ||
2150 * !locally_owned_water_dofs.is_element(dst));
2152 *
my_local_flow_routing.erase(it, my_local_flow_routing.end());
2156 * Pre-compute the write buffers that will be used in vmult(). These
2157 * buffers store the source and destination indices
for efficient batch
2158 * operations on the vector. This approach is more efficient than looping
2159 * over the flow routing pairs in each vmult() call because it allows us
2160 * to look up many vector entries all at once, rather than having to
2161 * translate between global and process-local indices
for each vector
2162 * entry we care about individually. The buffers are computed once in the
2163 * constructor and reused
for every matrix-vector product call.
2166 *
write_buffer_source_indices.resize(my_local_flow_routing.
size());
2167 *
write_buffer_indices.resize(my_local_flow_routing.
size());
2168 *
write_buffer_values.resize(my_local_flow_routing.
size());
2169 *
unsigned int index = 0;
2170 *
for (const auto &[src, dst] : my_local_flow_routing)
2172 *
write_buffer_source_indices[
index] = src;
2173 *
write_buffer_indices[
index] = dst;
2179 *
vmult(
typename VectorType::BlockType &y,
2180 *
const typename VectorType::BlockType &x)
const
2182 *
x_with_ghosts = x;
2186 * The src->dst relationship defines the
matrix via an entry
2187 * of +1 in the (src,src) position, and a -1 in the
2188 * (dst,src) position -- i.e., each entry in the src->dst
2189 * map defines a column of the
matrix.
2194 * that by setting y=I*x
2201 * Then we need to add to the locally-owned elements of the y vector
2202 * by multiplying the x vector with the -1's of
matrix.
2203 * This means that we need to
loop over all elements of the
2204 * map and determine whether the row
value of the entries
2205 * mentioned above are in the locally owned range:
2209 * Rather than looping over each (src, dst) pair and updating y(dst)
2210 * individually, we use a more efficient
vectorized approach via write
2212 * extract_subvector_to(), negate them, and then add them to the
2213 * destination vector
using array-based operations. This is much faster
2214 * than
scalar operations because it allows better use of the CPU
's
2215 * vectorization capabilities. The equivalent but slower code would be:
2216 * for (const auto &[src, dst] : my_local_flow_routing)
2217 * y(dst) -= x_with_ghosts(src);
2220 * x_with_ghosts.extract_subvector_to(write_buffer_source_indices,
2221 * write_buffer_values);
2222 * for (auto &v : write_buffer_values)
2224 * y.add(write_buffer_indices, write_buffer_values);
2226 * y.compress(VectorOperation::add);
2230 * mutable typename VectorType::BlockType x_with_ghosts;
2231 * mutable std::vector<types::global_dof_index> write_buffer_source_indices;
2232 * mutable std::vector<types::global_dof_index> write_buffer_indices;
2233 * mutable std::vector<PetscScalar> write_buffer_values;
2235 * std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
2236 * my_local_flow_routing;
2243 * <a name="parallel_flow_routing.cc-ParallelFlowRouterFlowRoutingPreconditioner"></a>
2244 * <h3>ParallelFlowRouter::FlowRoutingPreconditioner</h3>
2248 * This class implements an efficient preconditioner for the matrix A=(I-F).
2249 * Since the matrix is triangular (after reordering DoFs from high to low
2250 * elevation), the preconditioner uses a triangular solve to approximate
2251 * the inverse of (I-F).
2255 * The triangular solve works by processing the DoFs in order. For each
2256 * equation i (corresponding to water flow at node i), we compute:
2257 * y_i = (x_i + y_i) / a_ii
2258 * where a_ii = 1 (the diagonal entries of (I-F) are all 1), and y_i
2259 * accumulates contributions from upstream nodes. We then update downstream
2260 * equations by adding y_i to their y values.
2264 * Because the matrix is triangular and we process nodes in order from high
2265 * to low elevation, each node's solution depends only on upstream (higher)
2266 * nodes, which have already been processed. This makes the triangular solve
2267 * very efficient and cache-friendly.
2271 * See the detailed comments in the vmult() function below
for a complete
2272 * explanation of how the triangular solve is implemented.
2275 *
class ParallelFlowRouter::FlowRoutingPreconditioner
2278 *
FlowRoutingPreconditioner(
2279 *
const IndexSet &locally_owned_water_dofs,
2280 *
const IndexSet &locally_relevant_water_dofs,
2281 *
const unsigned int water_dofs_offset,
2284 *
: my_local_flow_routing(local_flow_routing)
2288 * We got the map from DoFs to downhill neighbors in global DoF
2289 * indices, but we need them in indices relative to the
second
2290 * vector block (or the (1,1) matrix block). So
shift, unless
2291 * the destination DoF is -1, indicating that
this source
2292 * DoF has no outlet (because it
's a depression in the DEM,
2293 * or because it's at the boundary).
2296 *
for (
auto &[src, dst] : my_local_flow_routing)
2298 *
src -= water_dofs_offset;
2299 *
Assert(locally_relevant_water_dofs.is_element(src),
2300 *
ExcInternalError());
2304 *
dst -= water_dofs_offset;
2305 *
Assert(locally_relevant_water_dofs.is_element(dst),
2306 *
ExcInternalError());
2312 * Unlike the
matrix itself, the preconditioner only looks at
2313 * the
diagonal blocks. Recall that
for each local routing src->dst,
2314 * we have entries in the (src,src) and (dst,src) position. Both
2315 * of these are in the same column,
'src'. One of the two entries
2320 * This means that
for a local routing to affect the
diagonal block,
2321 * we have to have
'src' be locally owned. That
's enough: if so,
2322 * at least the diagonal entry (and perhaps also the other one) is
2323 * in the locally owned diagonal block of the matrix.
2327 * So erase all others so that we don't have to
check this during
2328 * the vmult() operation:
2331 *
Assert(locally_owned_water_dofs.is_contiguous(), ExcInternalError());
2333 *
std::remove_if(my_local_flow_routing.
begin(),
2334 *
my_local_flow_routing.
end(),
2335 *
[&locally_owned_water_dofs](
2336 *
const
std::pair<
types::global_dof_index,
2337 *
types::global_dof_index> &src_dst) {
2339 *
return !locally_owned_water_dofs.is_element(src);
2341 *
my_local_flow_routing.erase(it, my_local_flow_routing.end());
2345 * At
this point, we should have
one routing
for each locally owned
2346 * DoF. Check that the number is right:
2350 *
locally_owned_water_dofs.n_elements());
2354 * One last step:
if we have a src->dst pair where we have already
2355 * made sure that
'src' is locally owned, then we know that the
2356 * (src,src) entry that results is in the locally owned diagonal
2357 * block. The
second matrix entry is (dst,src), which may or may
2358 * not be in that
diagonal block, depending on whether
'dst' is
2360 *
if 'src' simply has no
downstream neighbor.
2361 * In the
first case, the (dst,src)
matrix entry is of no concern to
2363 * In other words, in neither
case is there an entry (dst,src) that
2364 * we need to deal with.
2368 * To make our work easier, we turn the
first into the
second case
2369 * so that in the vmult() function we need not test
for inclusion
2370 * of 'dst' in the index
set of locally owned DoFs, but just compare
2371 * with
numbers::invalid_dof_index.
2374 *
for (auto &src_dst : my_local_flow_routing)
2376 *
(locally_owned_water_dofs.is_element(src_dst.
second) == false))
2448 *
vmult(typename VectorType::BlockType &y,
2449 *
const typename VectorType::BlockType &x) const
2452 *
for (
const auto &[src, dst] : my_local_flow_routing)
2456 * Solve the
'src'th equation. In the notation from
2457 * above, this reads as
2459 * which with this function
's variable names translates to
2460 * the following, storing the result in a temporary variable
2464 * const double yk = (y(src) += x(src));
2468 * Update a downstream entry if necessary. Again, in the notation
2469 * from above, this reads as
2471 * and so is the following:
2474 * if (dst != numbers::invalid_dof_index)
2477 * y.compress(VectorOperation::add);
2481 * std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
2482 * my_local_flow_routing;
2489 * <a name="parallel_flow_routing.cc-ParallelFlowRouterIplusminusXMatrixBase"></a>
2490 * <h3>ParallelFlowRouter::IplusminusXMatrixBase</h3>
2494 * This is a base class for matrix-free operators representing matrices of the
2495 * form (I ± X), where X is related to the flow routing matrix.
2499 * The class uses a block partitioning of the matrix and DoFs:
2500 * - Diagonal block (locally owned rows and columns)
2501 * - R matrix (locally owned rows, but columns on other processes)
2505 * The vmult() function computes y=(I ± X)*x by splitting the computation into
2506 * contributions from the locally owned part and the off-process part.
2510 * The derived classes IplusXMatrix and IminusXMatrix implement the two
2511 * variants of the actual operator.
2514 * class ParallelFlowRouter::IplusminusXMatrixBase
2517 * IplusminusXMatrixBase(
2518 * const IndexSet &locally_owned_water_dofs,
2519 * const IndexSet &locally_relevant_water_dofs,
2520 * const MPI_Comm mpi_communicator,
2521 * const unsigned int water_dofs_offset,
2522 * const std::vector<std::pair<types::global_dof_index,
2523 * types::global_dof_index>> &local_flow_routing,
2524 * const FlowRoutingPreconditioner &flow_routing_preconditioner)
2525 * : tmp(locally_owned_water_dofs, mpi_communicator)
2526 * , x_with_ghosts(locally_owned_water_dofs,
2527 * locally_relevant_water_dofs,
2529 * , my_local_flow_routing(local_flow_routing)
2530 * , flow_routing_preconditioner(flow_routing_preconditioner)
2534 * We got the map from DoFs to downhill neighbors in global DoF
2535 * indices, but we need them in indices relative to the second
2536 * vector block (or the (1,1) matrix block). So shift, unless
2537 * the destination DoF is -1, indicating that this source
2538 * DoF has no outlet (because it's a depression in the DEM,
2539 * or because it
's at the boundary).
2542 * for (auto &[src, dst] : my_local_flow_routing)
2544 * src -= water_dofs_offset;
2545 * Assert(locally_relevant_water_dofs.is_element(src),
2546 * ExcInternalError());
2548 * if (dst != numbers::invalid_dof_index)
2550 * dst -= water_dofs_offset;
2551 * Assert(locally_relevant_water_dofs.is_element(dst),
2552 * ExcInternalError());
2558 * Compared to the vmult() function of the FlowRoutingMatrix, where
2559 * we needed all entries that are in the locally owned rows of the
2560 * matrix, for the current matrix all we need are those entries that
2561 * are in the locally owned rows *but not in the locally owned columns*.
2562 * As a consequence, delete not only everything that's not in locally
2563 * owned rows, but *also* those *are* in locally owned columns:
2566 *
Assert(locally_owned_water_dofs.is_contiguous(), ExcInternalError());
2568 *
std::remove_if(my_local_flow_routing.begin(),
2569 *
my_local_flow_routing.end(),
2570 *
[&locally_owned_water_dofs](
2573 * const types::global_dof_index from = src_dst.first;
2574 * const types::global_dof_index to = src_dst.second;
2575 * return (locally_owned_water_dofs.is_element(from) ||
2576 * ((to == numbers::invalid_dof_index) ||
2577 * !locally_owned_water_dofs.is_element(to)));
2579 *
my_local_flow_routing.erase(it, my_local_flow_routing.end());
2584 * block that spans the whole
matrix, and so there is literally
nothing
2585 * left
for B. Make sure that is in fact
true.
2589 *
(my_local_flow_routing.size() == 0),
2590 *
ExcInternalError());
2594 * Pre-compute the write buffers
for efficient batch vector operations
2595 * in the vmult() function. This approach avoids repeated allocations and
2596 * allows vectorized extraction and addition operations on the flow
2597 * routing pairs, which is significantly faster than processing them
2601 *
write_buffer_source_indices.resize(my_local_flow_routing.
size());
2602 *
write_buffer_indices.resize(my_local_flow_routing.
size());
2603 *
write_buffer_values.resize(my_local_flow_routing.
size());
2604 *
unsigned int index = 0;
2605 *
for (const auto &[src, dst] : my_local_flow_routing)
2607 *
write_buffer_source_indices[
index] = src;
2608 *
write_buffer_indices[
index] = dst;
2614 *
mutable typename VectorType::BlockType tmp;
2615 *
mutable typename VectorType::BlockType x_with_ghosts;
2616 *
mutable std::vector<types::global_dof_index> write_buffer_source_indices;
2617 *
mutable std::vector<types::global_dof_index> write_buffer_indices;
2618 *
mutable std::vector<PetscScalar> write_buffer_values;
2620 *
std::vector<std::pair<types::global_dof_index, types::global_dof_index>>
2621 *
my_local_flow_routing;
2622 *
const FlowRoutingPreconditioner &flow_routing_preconditioner;
2629 * <a name=
"parallel_flow_routing.cc-ParallelFlowRouterIplusXMatrix"></a>
2630 * <h3>ParallelFlowRouter::IplusXMatrix</h3>
2634 * This
class represents the
matrix (I + X), where X is derived from the
2635 * flow routing
matrix. It is used in certain implicit time-stepping schemes.
2636 * The
class simply inherits from IplusminusXMatrixBase.
2639 *
class ParallelFlowRouter::IplusXMatrix
2640 *
:
public ParallelFlowRouter::IplusminusXMatrixBase
2643 *
IplusXMatrix(
const IndexSet &locally_owned_water_dofs,
2644 *
const IndexSet &locally_relevant_water_dofs,
2646 *
const unsigned int water_dofs_offset,
2647 *
const std::vector<
2648 *
std::pair<types::global_dof_index, types::global_dof_index>>
2649 *
&local_flow_routing,
2650 *
const FlowRoutingPreconditioner &flow_routing_preconditioner)
2651 *
: IplusminusXMatrixBase(locally_owned_water_dofs,
2652 *
locally_relevant_water_dofs,
2654 *
water_dofs_offset,
2655 *
local_flow_routing,
2656 *
flow_routing_preconditioner)
2661 *
vmult(
typename VectorType::BlockType &y,
2662 *
const typename VectorType::BlockType &x)
const
2666 * Start by importing ghost entries:
2669 *
x_with_ghosts = x;
2673 * The from->to relationship defines the
matrix via an entry
2674 * of +1 in the (from,from) position, and a -1 in the
2675 * (to,from) position -- i.e., each entry in the from->to
2676 * map defines a column of the
matrix.
2680 * The +1s all lie outside the D matrix, so we need not think
2681 * about these entries at all, and we can start with a zero vector:
2688 * Then we need to add to the locally-owned elements of the dst vector
2689 * by multiplying the src vector with the -1's of
matrix.
2690 * This means that we need to
loop over all elements of the
2691 * map and determine whether the row
value of the entries
2692 * mentioned above are in the locally owned range:
2696 * We use the write buffer optimization to efficiently compute the
2699 *
const std::pair<types::global_dof_index, types::global_dof_index>
2700 * locally_owned_range = tmp.local_range();
2701 *
for (
const auto &[src, dst] : my_local_flow_routing)
2704 * (dst >= locally_owned_range.first) &&
2705 * (dst < locally_owned_range.second),
2706 * ExcInternalError());
2707 * tmp(dst) -= x_with_ghosts(src);
2709 * However,
this version is much more efficient because it uses array
2711 * operations in a
loop.
2714 *
x_with_ghosts.extract_subvector_to(write_buffer_source_indices,
2715 *
write_buffer_values);
2716 *
for (
auto &v : write_buffer_values)
2718 *
tmp.add(write_buffer_indices, write_buffer_values);
2724 * We have now computed tmp=R*x. Let
's apply D^{-1} to it
2728 * flow_routing_preconditioner.vmult(y, tmp);
2732 * Finally, we need to add to it src so that we get (I+X)*src:
2743 * <a name="parallel_flow_routing.cc-ParallelFlowRouterIminusXMatrix"></a>
2744 * <h3>ParallelFlowRouter::IminusXMatrix</h3>
2748 * This class represents the matrix (I - X), where X is derived from the
2749 * flow routing matrix. The implementation is very similar to IplusXMatrix,
2750 * but with the opposite sign for X.
2753 * class ParallelFlowRouter::IminusXMatrix
2754 * : public ParallelFlowRouter::IplusminusXMatrixBase
2757 * IminusXMatrix(const IndexSet &locally_owned_water_dofs,
2758 * const IndexSet &locally_relevant_water_dofs,
2759 * const MPI_Comm mpi_communicator,
2760 * const unsigned int water_dofs_offset,
2761 * const std::vector<
2762 * std::pair<types::global_dof_index, types::global_dof_index>>
2763 * &local_flow_routing,
2764 * const FlowRoutingPreconditioner &flow_routing_preconditioner)
2765 * : IplusminusXMatrixBase(locally_owned_water_dofs,
2766 * locally_relevant_water_dofs,
2768 * water_dofs_offset,
2769 * local_flow_routing,
2770 * flow_routing_preconditioner)
2775 * vmult(typename VectorType::BlockType &y,
2776 * const typename VectorType::BlockType &x) const
2778 * x_with_ghosts = x;
2781 * x_with_ghosts.extract_subvector_to(write_buffer_source_indices,
2782 * write_buffer_values);
2783 * tmp.add(write_buffer_indices, write_buffer_values);
2785 * tmp.compress(VectorOperation::add);
2789 * We have now computed tmp=-R*x. Let's
apply D^{-1} to it
2793 *
flow_routing_preconditioner.vmult(y, tmp);
2797 * Finally, we need to add to it src so that we get (I-X)*src:
2808 * <a name=
"parallel_flow_routing.cc-ParallelFlowRouterassemble_matrix_free_operators"></a>
2809 * <h3>ParallelFlowRouter::assemble_matrix_free_operators()</h3>
2813 * Rather than assembling the full
matrix (which could be very large in
2814 *
parallel),
this function creates
"matrix-free" operators that compute
2815 *
matrix-vector products implicitly. We will then be able to compare
2816 * all of these approaches.
2820 * Specifically,
this function creates three key objects:
2821 * 1.
A FlowRoutingMatrix that represents the
matrix F described in
2822 * assemble_system() above.
2823 * 2. A FlowRoutingPreconditioner that approximates the inverse of (I - F)
2824 * using a fast triangular solve.
2825 * 3. A pair of matrices I+X and I-X.
2829 * These matrix-free operators are used in the solve() function to perform
2830 * iterative linear solves without ever explicitly storing the full matrix.
2834 *
ParallelFlowRouter::assemble_matrix_free_operators()
2837 *
"Solver 2: Assemble matrix-free operators");
2838 *
pcout <<
"Assembling matrix-free operators... " << std::endl;
2840 *
flow_routing_matrix = std::make_unique<const FlowRoutingMatrix>(
2841 *
locally_owned_partitioning[1],
2842 *
locally_relevant_partitioning[1],
2844 *
locally_relevant_partitioning[0].
size(),
2845 *
local_flow_routing);
2847 *
flow_routing_preconditioner =
2848 *
std::make_unique<const FlowRoutingPreconditioner>(
2849 *
locally_owned_partitioning[1],
2850 *
locally_relevant_partitioning[1],
2851 *
locally_relevant_partitioning[0].
size(),
2852 *
local_flow_routing);
2854 *
I_plus_X_matrix = std::make_unique<const IplusXMatrix>(
2855 *
locally_owned_partitioning[1],
2856 *
locally_relevant_partitioning[1],
2858 *
locally_relevant_partitioning[0].
size(),
2859 *
local_flow_routing,
2860 *
*flow_routing_preconditioner);
2862 *
I_minus_X_matrix = std::make_unique<const IminusXMatrix>(
2863 *
locally_owned_partitioning[1],
2864 *
locally_relevant_partitioning[1],
2866 *
locally_relevant_partitioning[0].
size(),
2867 *
local_flow_routing,
2868 *
*flow_routing_preconditioner);
2876 * <a name=
"parallel_flow_routing.cc-ParallelFlowRoutersolve"></a>
2877 * <h3>ParallelFlowRouter::solve()</h3>
2881 * This function solves the linear system assembled in assemble_system() to
2882 * find the steady-state water flow rates at all nodes on the mesh. Since we
2883 * have renumbered the degrees of freedom from high to low elevation, the
2884 * system matrix is lower triangular, and we can solve it very efficiently.
2888 * The function uses an iterative Richardson solver with a preconditioner
2889 * derived from the triangular structure of the matrix. Because the matrix
2890 * is triangular, a simple SOR preconditioner with a relaxation factor of
2891 * 1.0 is equivalent to the triangular solve.
2895 * The function then solves the system three more times using the matrix-free
2896 * operators defined above, and compares the results to verify that they all
2897 * give the same solution.
2901 *
ParallelFlowRouter::solve()
2903 *
pcout <<
"Solving for global water routing... " << std::endl;
2907 * ---------------- Solve
matrix-based -------------------------
2910 *
VectorType completely_distributed_solution_matrix_based(
2911 *
locally_owned_partitioning, mpi_communicator);
2914 *
"Solver 1: Solve for water matrix-based");
2917 *
1e-6 * system_rhs.block(1).l2_norm());
2921 *
preconditioner.
initialize(system_matrix.block(1, 1));
2923 *
solver.solve(system_matrix.block(1, 1),
2924 *
completely_distributed_solution_matrix_based.block(1),
2925 *
system_rhs.block(1),
2928 *
pcout <<
" Solved matrix-based in " << solver_control.last_step()
2929 *
<<
" iterations." << std::endl;
2935 * ---------------- Solve
matrix-free -------------------------
2938 *
VectorType completely_distributed_solution_matrix_free(
2939 *
locally_owned_partitioning, mpi_communicator);
2942 *
"Solver 2: Solve for water matrix-free");
2945 *
1e-6 * system_rhs.block(1).l2_norm());
2948 *
solver.solve(*flow_routing_matrix,
2949 *
completely_distributed_solution_matrix_free.block(1),
2950 *
system_rhs.block(1),
2951 *
*flow_routing_preconditioner);
2953 *
pcout <<
" Solved matrix-free in " << solver_control.last_step()
2954 *
<<
" iterations." << std::endl;
2959 * ---------------- Solve via (I+X)x=D^{-1}
b -------------------------
2962 *
VectorType completely_distributed_solution_IplusX(
2963 *
locally_owned_partitioning, mpi_communicator);
2967 *
typename VectorType::BlockType Dinv_times_rhs(
2968 *
locally_owned_partitioning[1], mpi_communicator);
2969 *
flow_routing_preconditioner->vmult(Dinv_times_rhs, system_rhs.block(1));
2972 *
1e-6 * Dinv_times_rhs.l2_norm());
2975 *
solver.solve(*I_plus_X_matrix,
2976 *
completely_distributed_solution_IplusX.block(1),
2980 *
pcout <<
" Solved I+X-based in " << solver_control.last_step()
2981 *
<<
" iterations." << std::endl;
2986 * ---------------- Solve via (I-X)(I+X)x=(I-X)D^{-1}
b -------------------
2989 *
VectorType completely_distributed_solution_IminusX_IplusX(
2990 *
locally_owned_partitioning, mpi_communicator);
2993 *
"Solver 4: Solve for water (I-X)(I+X)");
2995 *
typename VectorType::BlockType Dinv_times_rhs(
2996 *
locally_owned_partitioning[1], mpi_communicator);
2997 *
flow_routing_preconditioner->vmult(Dinv_times_rhs, system_rhs.block(1));
3000 *
1e-6 * Dinv_times_rhs.l2_norm());
3003 *
solver.solve(*I_plus_X_matrix,
3004 *
completely_distributed_solution_IminusX_IplusX.block(1),
3006 *
*I_minus_X_matrix);
3008 *
pcout <<
" Solved (I-X)(I+X)-based in " << solver_control.last_step()
3009 *
<<
" iterations." << std::endl;
3012 *
locally_relevant_solution.block(1) =
3013 *
completely_distributed_solution_matrix_free.block(1);
3017 * ----------- Now make sure the solutions agree:
3020 *
completely_distributed_solution_matrix_free -=
3021 *
completely_distributed_solution_matrix_based;
3022 *
pcout <<
" Relative error between matrix-based and matrix-free: "
3023 *
<< completely_distributed_solution_matrix_free.l2_norm() /
3024 *
completely_distributed_solution_matrix_based.l2_norm()
3027 *
completely_distributed_solution_IplusX -=
3028 *
completely_distributed_solution_matrix_based;
3029 *
pcout <<
" Relative error between matrix-based and I+X solution: "
3030 *
<< completely_distributed_solution_IplusX.l2_norm() /
3031 *
completely_distributed_solution_matrix_based.l2_norm()
3034 *
completely_distributed_solution_IminusX_IplusX -=
3035 *
completely_distributed_solution_matrix_based;
3036 *
pcout <<
" Relative error between matrix-based and (I-X)(I+X) solution: "
3037 *
<< completely_distributed_solution_IminusX_IplusX.l2_norm() /
3038 *
completely_distributed_solution_matrix_based.l2_norm()
3046 * <a name=
"parallel_flow_routing.cc-ParallelFlowRoutercheck_conservation_for_waterflow_system"></a>
3047 * <h3>ParallelFlowRouter::check_conservation_for_waterflow_system()</h3>
3051 * This function verifies that the solution satisfies the principle of mass
3052 * conservation
for water. Specifically, it checks that the total water input
3053 * (from rainfall) equals the total water output (flowing out of the domain).
3057 * The function computes:
3058 * - The total water input by integrating the rainfall rate over the domain
3059 * - The total water output by summing the water flow rates at boundary nodes
3063 * If these two quantities differ by more than a small tolerance (currently
3064 * 1%), the function
throws an exception, indicating a problem with the
3069 * This is a valuable diagnostic
check: in a correctly formulated and solved
3070 * system, mass should be strictly conserved (up to numerical errors). If
3071 * conservation is violated, it indicates an error in problem setup, assembly,
3076 *
ParallelFlowRouter::check_conservation_for_waterflow_system(
3077 *
const VectorType &solution)
3081 *
const QGauss<dim> quadrature_formula(fe.degree + 1);
3083 *
quadrature_formula,
3087 *
const RainFallRate<spacedim> rainfall_rate;
3088 *
double input_from_rain_rate = 0.0;
3089 *
for (
const auto &cell : dof_handler.active_cell_iterators())
3090 *
if (cell->is_locally_owned())
3092 *
fe_values.
reinit(cell);
3094 *
for (
const unsigned int q : fe_values.quadrature_point_indices())
3095 *
input_from_rain_rate +=
3096 *
rainfall_rate.
value(fe_values.quadrature_point(q)) *
3100 *
input_from_rain_rate =
3106 * Now also
check the outflow. Water flows out of the domain at points
3107 * that (i) are at the boundary, and (ii) have no lower neighbors. We
3108 * have built the DEM so that it has no local depressions in the
3109 * interior of the domain, so we only have to
check (ii). In the
3110 * local water routing table,
this is indicated by (src->dst) pairs
3111 * where 'dst' is an invalid_dof_index. The only thing we have to pay
3112 * attention to is that we only count the locally owned DoFs to
3113 * avoid
double-counting:
3116 *
double water_outflow_rate = 0;
3117 *
for (const auto &[src, dst] : local_flow_routing)
3118 *
if ((dst ==
numbers::invalid_dof_index) &&
3119 *
dof_handler.locally_owned_dofs().is_element(src))
3120 *
water_outflow_rate += solution(src);
3122 *
water_outflow_rate =
3123 *
Utilities::
MPI::sum(water_outflow_rate, mpi_communicator);
3125 *
const
double error_abs =
3126 *
std::abs(input_from_rain_rate - water_outflow_rate);
3127 *
const
double error_rel =
std::abs(error_abs / input_from_rain_rate);
3128 *
pcout << "Conservation check (water)" <<
std::endl
3129 *
<< " Input: " << input_from_rain_rate <<
std::endl
3130 *
<< " Output: " << water_outflow_rate <<
std::endl
3131 *
<< " Relative error: " << error_rel <<
std::endl;
3134 *
ExcMessage("Conservation of water rate not satisfied."));
3141 * <a name="parallel_flow_routing.cc-ParallelFlowRouteroutput_results"></a>
3142 * <h3>ParallelFlowRouter::output_results()</h3>
3146 * This function writes the solution to VTU output files
for visualization.
3147 * The function also adds the subdomain ID of each cell, which is useful
for
3148 * visualizing the
parallel partitioning of the mesh across
MPI processes.
3149 * This can help verify that the partitioning is reasonably balanced.
3153 * In
parallel, each process writes its local portion of the solution to a
3154 * separate file, and a master PVTU file is created that ties all the pieces
3155 * together
for visualization in tools like ParaView or VisIt.
3159 *
ParallelFlowRouter::output_results()
3162 *
pcout <<
"Writing output... " << std::flush;
3164 *
const std::vector<std::string> solution_names = {
"elevation",
3165 *
"water_flow_rate"};
3166 *
const std::vector<DataComponentInterpretation::DataComponentInterpretation>
3167 *
data_component_interpretation = {
3173 *
data_out.add_data_vector(locally_relevant_solution,
3176 *
data_component_interpretation);
3179 *
for (
unsigned int i = 0; i < subdomain.size(); ++i)
3180 *
subdomain(i) = triangulation.locally_owned_subdomain();
3181 *
data_out.add_data_vector(subdomain,
"subdomain");
3183 *
data_out.build_patches();
3185 *
data_out.write_vtu_with_pvtu_record(
3186 *
"./",
"solution", 0, mpi_communicator, 2);
3194 * <a name=
"parallel_flow_routing.cc-ParallelFlowRouterrun"></a>
3195 * <h3>ParallelFlowRouter::run()</h3>
3199 * This is the
main entry
point for the flow routing solver. It orchestrates
3200 * all the steps needed to solve the water flow routing problem:
3204 * 1. Create the computational mesh (`make_grid()`)
3205 * 2. Set up the finite element spaces and DoF numbering (`setup_dofs()`)
3206 * 3. Interpolate the digital elevation model onto the mesh
3207 * (`interpolate_initial_elevation()`)
3208 * 4. Renumber DoFs so water flows from high to low elevation
3209 * (`sort_dofs_high_to_low()`)
3210 * 5. Determine which downhill neighbor each node flows to
3211 * (`compute_local_flow_routing()`)
3212 * 6. Assemble the linear system
for steady-state flow (`assemble_system()`)
3213 * 7. Set up matrix-free operators
for efficient solving
3214 * (`assemble_matrix_free_operators()`)
3215 * 8. Solve the linear system (`solve()`)
3216 * 9. Check conservation of water mass
3217 * (`check_conservation_for_waterflow_system()`)
3218 * 10. Write output
for visualization (`output_results()`)
3219 * 11. Print performance statistics (`computing_timer.print_summary()`)
3223 * The sequence of these steps reflects the logical flow of the algorithm:
3224 * we
first set up the geometry and DoFs, then compute the flow routing
3225 * connectivity, then
assemble and solve the linear system, and finally
3226 * perform validation and output.
3230 *
ParallelFlowRouter::
run()
3234 *
interpolate_initial_elevation();
3235 *
sort_dofs_high_to_low();
3237 *
compute_local_flow_routing();
3239 *
assemble_system();
3240 *
assemble_matrix_free_operators();
3244 *
check_conservation_for_waterflow_system(locally_relevant_solution);
3245 *
if (generate_graphical_output)
3250 * Print the time taken
for each section of the code,
first in the summary
3251 * table and then as individual
numbers in
one line. This is useful
for
3252 * creating graphs of
run times.
3255 *
computing_timer.print_summary();
3256 *
pcout <<
"Times per section: ";
3257 *
for (
const auto &[name, time] :
3258 *
computing_timer.get_summary_data(
TimerOutput::total_wall_time))
3259 *
pcout << time <<
' ';
3260 *
pcout << std::endl;
3269 * <a name=
"parallel_flow_routing.cc-Themainfunction"></a>
3270 * <h3>The
main() function</h3>
3274 * The
main function of the program is quite simple. It initializes
MPI for
3275 *
parallel execution, creates a ParallelFlowRouter
object that registers all
3276 * necessary parameters, and then runs the flow routing algorithm.
3280 * The parameters can be provided in a file (by passing the filename as a
3281 * command-line argument) or are left at their defaults.
3285 * The function includes basic error handling: if an exception occurs during
3286 * the computation, it prints an error message and exits gracefully with a
3287 * non-zero return code.
3291 *
main(
int argc,
char *argv[])
3295 *
using namespace dealii;
3301 * Create the problem object (
this registers parameters with
3305 *
ParallelFlowRouting::ParallelFlowRouter problem;
3309 * Parse parameters from file
if provided
3319 *
catch (std::exception &exc)
3321 *
std::cerr << std::endl
3323 *
<<
"----------------------------------------------------"
3325 *
std::cerr <<
"Exception on processing: " << std::endl
3326 *
<< exc.what() << std::endl
3327 *
<<
"Aborting!" << std::endl
3328 *
<<
"----------------------------------------------------"
3335 *
std::cerr << std::endl
3337 *
<<
"----------------------------------------------------"
3339 *
std::cerr <<
"Unknown exception!" << std::endl
3340 *
<<
"Aborting!" << std::endl
3341 *
<<
"----------------------------------------------------"
* * for(const auto &cell :triangulation.active_cell_iterators())
* * int main(int argc, char **argv)
* x_component_mask set(0, true)
* * * struct InterferenceTaperTransform *
***mech_lbc_system increment_interpolation_handlers push_back(scale_z_handler)
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
IndexSet get_view(const size_type begin, const size_type end) const
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
static void initialize(const std::string &filename="", const std::string &output_filename="", const ParameterHandler::OutputStyle output_style_for_output_filename=ParameterHandler::Short, ParameterHandler &prm=ParameterAcceptor::prm, const ParameterHandler::OutputStyle output_style_for_filename=ParameterHandler::DefaultStyle)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
Expression atan2(const Expression &y, const Expression &x)
Expression operator>(const Expression &lhs, const Expression &rhs)
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void downstream(DoFHandler< dim, spacedim > &dof_handler, const Tensor< 1, spacedim > &direction, const bool dof_wise_renumbering=false)
void random(DoFHandler< dim, spacedim > &dof_handler)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell< dim > &reference_cell)
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
constexpr types::blas_int zero
constexpr types::blas_int one
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
* * * ScaleZFunction< dim, Number, components >::ScaleZFunction * component(component)
* * if(update_pressure &update_flags) * compute_pressure(constitutive_request
* * * * std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters const
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, ShapeDataMemorySpace > shape_data, const ViewTypeIn in, ViewTypeOut out)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
std::vector< T > all_gather(const MPI_Comm comm, const T &object_to_send)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
T reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
T broadcast(const MPI_Comm comm, const T &object_to_send, const unsigned int root_process=0)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
constexpr types::global_dof_index invalid_dof_index
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
1.9.8