Merge pull request #5603 from gassmoeller/split_heat_flux_output

Improve heat flux map postprocessor and dynamic topography postprocessor
geodynamics · Mar 30, 2024 · 946ecab · 946ecab
2 parents 6adb11c + db58e71
commit 946ecab
Show file tree

Hide file tree

Showing 27 changed files with 208 additions and 175 deletions.
diff --git a/doc/modules/changes/20240315_gassmoeller2 b/doc/modules/changes/20240315_gassmoeller2
@@ -0,0 +1,8 @@
+Changed: The postprocessor 'heat flux map' now outputs the heat flux values in
+the same format as the 'dynamic topography' postprocessor does for dynamic
+topography. This in particular means instead of writing all values into a
+single file 'heat flux map' now writes one file per boundary. Both
+postprocessors now take care to not include spaces in the data file
+headers so that separating the header columns by spaces works correctly.
+<br>
+(Rene Gassmoeller, 2024/03/15)
diff --git a/include/aspect/postprocess/dynamic_topography.h b/include/aspect/postprocess/dynamic_topography.h
@@ -87,9 +87,12 @@ namespace aspect
       private:
         /**
          * Output the dynamic topography solution to
-         * a file.
+         * a file for a given boundary id. All the values
+         * in @p position_and_topography are written to a file
+         * with a file name determined by @p boundary_id.
          */
-        void output_to_file(bool upper, std::vector<std::pair<Point<dim>, double>> &values);
+        void output_to_file(const types::boundary_id boundary_id,
+                            const std::vector<std::pair<Point<dim>, double>> &position_and_topography);
 
         /**
          * A vector which stores the surface stress values calculated

diff --git a/include/aspect/postprocess/heat_flux_map.h b/include/aspect/postprocess/heat_flux_map.h
@@ -88,6 +88,19 @@ namespace aspect
          */
         std::pair<std::string,std::string>
         execute (TableHandler &statistics) override;
+
+      private:
+        /**
+         * Output the heat flux density for the boundary determined
+         * by @p boundary_id to a file. The heat flux density is
+         * handed over in the vector @p heat_flux_and_area. This vector
+         * is expected to be of the structure described for the return value
+         * of the function compute_heat_flux_through_boundary_faces() and
+         * only the values at the faces of the given @p boundary_id are
+         * written to the file.
+         */
+        void output_to_file(const types::boundary_id boundary_id,
+                            const std::vector<std::vector<std::pair<double, double>>> &heat_flux_and_area);
     };
   }
 }

diff --git a/source/postprocess/dynamic_topography.cc b/source/postprocess/dynamic_topography.cc
@@ -365,9 +365,9 @@ namespace aspect
 
       // Possibly output the result to file.
       if (output_surface)
-        output_to_file(true, stored_values_surface);
+        output_to_file(top_boundary_id, stored_values_surface);
       if (output_bottom)
-        output_to_file(false, stored_values_bottom);
+        output_to_file(bottom_boundary_id, stored_values_bottom);
 
       return std::pair<std::string,std::string>("Computing dynamic topography", "");
     }
@@ -411,32 +411,35 @@ namespace aspect
      */
     template <int dim>
     void
-    DynamicTopography<dim>::output_to_file(const bool upper,
-                                           std::vector<std::pair<Point<dim>,
-                                           double>> &stored_values)
+    DynamicTopography<dim>::output_to_file(const types::boundary_id boundary_id,
+                                           const std::vector<std::pair<Point<dim>,double>> &position_and_topography)
     {
+      // get boundary name and avoid spaces for file output
+      std::string boundary_name = this->get_geometry_model().translate_id_to_symbol_name(boundary_id);
+      std::replace(boundary_name.begin(), boundary_name.end(), ' ', '_');
+
       std::ostringstream output;
 
       // On processor 0, write header lines
       if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) == 0)
         {
           output << "# "
                  << ((dim==2)? "x y " : "x y z ")
-                 << (upper ? "surface topography" : "bottom topography") << std::endl;
+                 << "dynamic_topography_" << boundary_name << std::endl;
         }
 
-      for (unsigned int i=0; i<stored_values.size(); ++i)
+      for (unsigned int i = 0; i < position_and_topography.size(); ++i)
         {
           output << std::setprecision(10)
-                 << stored_values[i].first
+                 << position_and_topography[i].first
                  << ' '
                  << std::setprecision(10)
-                 << stored_values[i].second
+                 << position_and_topography[i].second
                  << std::endl;
         }
 
       std::string filename = this->get_output_directory() +
-                             (upper ? "dynamic_topography_surface." : "dynamic_topography_bottom.") +
+                             "dynamic_topography_" + boundary_name + "." +
                              Utilities::int_to_string(this->get_timestep_number(), 5);
       if (this->get_parameters().run_postprocessors_on_nonlinear_iterations)
         filename.append("." + Utilities::int_to_string (this->get_nonlinear_iteration(), 4));

diff --git a/source/postprocess/heat_flux_map.cc b/source/postprocess/heat_flux_map.cc
@@ -460,30 +460,46 @@ namespace aspect
       std::vector<std::vector<std::pair<double, double>>> heat_flux_and_area =
         internal::compute_heat_flux_through_boundary_faces (*this);
 
-      // have a stream into which we write the data. the text stream is then
-      // later sent to processor 0
+      const auto boundary_ids = this->get_geometry_model().get_used_boundary_indicators();
+      for (const auto &boundary_id: boundary_ids)
+        if (this->get_geometry_model().translate_id_to_symbol_name(boundary_id) == "top" ||
+            this->get_geometry_model().translate_id_to_symbol_name(boundary_id) == "bottom")
+          output_to_file(boundary_id, heat_flux_and_area);
+
+      std::string placeholder_name = this->get_output_directory() + "heat_flux." + Utilities::int_to_string(this->get_timestep_number(), 5);
+      if (this->get_parameters().run_postprocessors_on_nonlinear_iterations)
+        placeholder_name.append("." + Utilities::int_to_string (this->get_nonlinear_iteration(), 4));
+
+      return std::pair<std::string,std::string>("Writing heat flux map",
+                                                placeholder_name);
+    }
+
+
+
+    template <int dim>
+    void
+    HeatFluxMap<dim>::output_to_file(const types::boundary_id boundary_id,
+                                     const std::vector<std::vector<std::pair<double, double>>> &heat_flux_and_area)
+    {
+      // get boundary name and avoid spaces for file output
+      std::string boundary_name = this->get_geometry_model().translate_id_to_symbol_name(boundary_id);
+      std::replace(boundary_name.begin(), boundary_name.end(), ' ', '_');
+
       std::ostringstream output;
 
       // write the file header. note that we only do so on processor 0
       if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) == 0)
         {
           output << "# "
                  << ((dim==2)? "x y" : "x y z")
-                 << " heat flux" << std::endl;
+                 << " heat_flux_" << boundary_name << std::endl;
         }
 
-      std::vector<std::pair<Point<dim>,double>> stored_values;
-
-      // loop over all of the surface cells and evaluate the heat flux
-      const types::boundary_id top_boundary_id = this->get_geometry_model().translate_symbolic_boundary_name_to_id("top");
-      const types::boundary_id bottom_boundary_id = this->get_geometry_model().translate_symbolic_boundary_name_to_id("bottom");
-
       for (const auto &cell : this->get_dof_handler().active_cell_iterators())
         if (cell->is_locally_owned())
-          for (const unsigned int f : cell->face_indices())
+          for (const unsigned int f: cell->face_indices())
             if (cell->at_boundary(f) &&
-                (cell->face(f)->boundary_id() == top_boundary_id ||
-                 cell->face(f)->boundary_id() == bottom_boundary_id))
+                (cell->face(f)->boundary_id() == boundary_id))
               {
                 // evaluate position of heat flow to write into output file
                 const bool respect_manifold = true;
@@ -492,33 +508,21 @@ namespace aspect
                 const double flux_density = heat_flux_and_area[cell->active_cell_index()][f].first /
                                             heat_flux_and_area[cell->active_cell_index()][f].second;
 
-                // store final position and heat flow
-                stored_values.emplace_back (midpoint_at_surface, flux_density);
+                output << std::setprecision(10)
+                       << midpoint_at_surface
+                       << ' '
+                       << std::setprecision(10)
+                       << flux_density
+                       << std::endl;
               }
 
-
-      // Write the solution to an output stream
-      for (unsigned int i=0; i<stored_values.size(); ++i)
-        {
-          output << stored_values[i].first
-                 << ' '
-                 << stored_values[i].second
-                 << std::endl;
-        }
-
-
-
       std::string filename = this->get_output_directory() +
-                             "heat_flux." +
+                             "heat_flux_" + boundary_name + "." +
                              Utilities::int_to_string(this->get_timestep_number(), 5);
       if (this->get_parameters().run_postprocessors_on_nonlinear_iterations)
         filename.append("." + Utilities::int_to_string (this->get_nonlinear_iteration(), 4));
 
-
       Utilities::collect_and_write_file_content(filename, output.str(), this->get_mpi_communicator());
-
-      return std::pair<std::string,std::string>("Writing heat flux map:",
-                                                filename);
     }
   }
 }

diff --git a/tests/dynamic_topography/dynamic_topography_bottom.00000 b/tests/dynamic_topography/dynamic_topography_bottom.00000
@@ -1,9 +1,9 @@
-# x y bottom topography
-112500 0 -20.35932195
-337500 0 -22.5880606
-562500 0 -13.12252722
-787500 0 55.68976695
-1012500 0 55.68976696
-1237500 0 -13.12252721
-1462500 0 -22.58806058
-1687500 0 -20.35932193
+# x y dynamic_topography_bottom
+112500 0 -20.35932198
+337500 0 -22.58806059
+562500 0 -13.12252715
+787500 0 55.68976686
+1012500 0 55.68976686
+1237500 0 -13.12252713
+1462500 0 -22.58806056
+1687500 0 -20.35932195
diff --git a/...pography/dynamic_topography_surface.00000 → ...c_topography/dynamic_topography_top.00000 b/...pography/dynamic_topography_surface.00000 → ...c_topography/dynamic_topography_top.00000
@@ -1,4 +1,4 @@
-# x y surface topography
+# x y dynamic_topography_top
 112500 700000 -37.56967558
 337500 700000 -37.24917518
 562500 700000 0.09274690416

diff --git a/tests/dynamic_topography2/dynamic_topography_bottom.00000 b/tests/dynamic_topography2/dynamic_topography_bottom.00000
@@ -1,4 +1,4 @@
-# x y bottom topography
+# x y dynamic_topography_bottom
 112500 0 -0.691340274
 337500 0 -0.4081226042
 562500 0 0.2240167169

diff --git a/...ography2/dynamic_topography_surface.00000 → ..._topography2/dynamic_topography_top.00000 b/...ography2/dynamic_topography_surface.00000 → ..._topography2/dynamic_topography_top.00000
@@ -1,4 +1,4 @@
-# x y surface topography
+# x y dynamic_topography_top
 112500 700000 -219.5206365
 337500 700000 -214.9051988
 562500 700000 -294.199505

diff --git a/tests/dynamic_topography3/dynamic_topography_bottom.00000 b/tests/dynamic_topography3/dynamic_topography_bottom.00000
@@ -1,4 +1,4 @@
-# x y z bottom topography
+# x y z dynamic_topography_bottom
 -1265819.559 -3241620.096 2.084595482e-11 -109.9679245
 -2324756.457 -2324756.457 1140795.701 -1130.427916
 -2324756.457 -2324756.457 -1140795.701 -1130.428354
@@ -28,7 +28,7 @@
 2324756.457 2324756.457 1140795.701 -1130.428763
 1265819.559 3241620.096 -1.518641958e-11 -109.9683258
 2324756.457 -1140795.701 -2324756.457 -922.8374217
-1265819.559 2.261096449e-11 -3241620.096 12.95300895
+1265819.559 2.261096449e-11 -3241620.096 12.95300896
 3241620.096 -2.335326054e-10 -1265819.559 172.706491
 2324756.457 1140795.701 -2324756.457 -922.8423924
 -1.131907049e-11 1265819.559 -3241620.096 -232.9154335

diff --git a/...ography3/dynamic_topography_surface.00000 → ..._topography3/dynamic_topography_top.00000 b/...ography3/dynamic_topography_surface.00000 → ..._topography3/dynamic_topography_top.00000
@@ -1,4 +1,4 @@
-# x y z surface topography
+# x y z dynamic_topography_top
 -2304664.576 -5901984.174 -1.717376213e-13 81.6990232
 -4232660.032 -4232660.032 2077034.932 -64.39460075
 -4232660.032 -4232660.032 -2077034.932 -64.39308593
@@ -38,7 +38,7 @@
 -4232660.032 4232660.032 -2077034.932 -64.39260648
 -5901984.174 2304664.576 -1.984763256e-10 81.69886784
 -2304664.576 5901984.174 -1.976176375e-10 81.69907751
--4232660.032 4232660.032 2077034.932 -64.3942961
+-4232660.032 4232660.032 2077034.932 -64.39429609
 -2304664.576 3.743880144e-11 -5901984.174 -158.7562861
 -4232660.032 -2077034.932 -4232660.032 -469.4021192
 -4232660.032 2077034.932 -4232660.032 -469.4152887

diff --git a/tests/dynamic_topography_back/dynamic_topography_bottom.00000 b/tests/dynamic_topography_back/dynamic_topography_bottom.00000
@@ -1,9 +1,9 @@
-# x y bottom topography
-112500 0 -20.35932195
-337500 0 -22.5880606
-562500 0 -13.12252722
-787500 0 55.68976695
-1012500 0 55.68976696
-1237500 0 -13.12252721
-1462500 0 -22.58806058
-1687500 0 -20.35932193
+# x y dynamic_topography_bottom
+112500 0 -20.35932198
+337500 0 -22.58806059
+562500 0 -13.12252715
+787500 0 55.68976686
+1012500 0 55.68976686
+1237500 0 -13.12252713
+1462500 0 -22.58806056
+1687500 0 -20.35932195
diff --git a/...pography/dynamic_topography_surface.00000 → ...ography_back/dynamic_topography_top.00000 b/...pography/dynamic_topography_surface.00000 → ...ography_back/dynamic_topography_top.00000
@@ -1,4 +1,4 @@
-# x y surface topography
+# x y dynamic_topography_top
 112500 700000 -37.56967558
 337500 700000 -37.24917518
 562500 700000 0.09274690416

diff --git a/tests/dynamic_topography_vs_surface_dynamic_topography/dynamic_topography_bottom.00000 b/tests/dynamic_topography_vs_surface_dynamic_topography/dynamic_topography_bottom.00000
@@ -1,9 +1,9 @@
-# x y bottom topography
-112500 0 -20.35932195
-337500 0 -22.5880606
-562500 0 -13.12252722
-787500 0 55.68976695
-1012500 0 55.68976696
-1237500 0 -13.12252721
-1462500 0 -22.58806058
-1687500 0 -20.35932193
+# x y dynamic_topography_bottom
+112500 0 -20.35932198
+337500 0 -22.58806059
+562500 0 -13.12252715
+787500 0 55.68976686
+1012500 0 55.68976686
+1237500 0 -13.12252713
+1462500 0 -22.58806056
+1687500 0 -20.35932195
diff --git a/...phy_back/dynamic_topography_surface.00000 → ...c_topography/dynamic_topography_top.00000 b/...phy_back/dynamic_topography_surface.00000 → ...c_topography/dynamic_topography_top.00000
@@ -1,4 +1,4 @@
-# x y surface topography
+# x y dynamic_topography_top
 112500 700000 -37.56967558
 337500 700000 -37.24917518
 562500 700000 0.09274690416

diff --git a/tests/geoid/dynamic_topography_bottom.00000 b/tests/geoid/dynamic_topography_bottom.00000
@@ -1,4 +1,4 @@
-# x y z bottom topography
+# x y z dynamic_topography_bottom
 -35498.28088 -196824.4702 -2.76471554e-13 -64957.61326
 -70785.07178 -184010.7447 33608.32387 -57553.88554
 -70785.07178 -184010.7447 -33608.32387 -57553.88558
@@ -21,10 +21,10 @@
 -139384.3779 -33822.9268 139384.3779 17280.58779
 -70785.07178 -33608.32387 184010.7447 46185.62438
 -35498.28088 2.394570686e-12 196824.4702 53281.72759
--109044.9905 8.994181261e-12 167657.9555 36339.18081
--70785.07178 33608.32387 184010.7447 46185.62439
+-109044.9905 1.06639181e-11 167657.9555 36339.18081
+-70785.07178 33608.32387 184010.7447 46185.62438
 -184010.7447 -33608.32387 70785.07178 -34581.30002
--167657.9555 2.303524681e-11 109044.9905 -5575.059541
+-167657.9555 2.422576705e-11 109044.9905 -5575.059541
 -196824.4702 -6.462091252e-12 35498.28088 -55980.07258
 -184010.7447 33608.32387 70785.07178 -34581.30001
 -139384.3779 33822.9268 139384.3779 17280.58779
@@ -34,7 +34,7 @@
 1.139428721e-12 -196824.4702 35498.28088 -55980.07254
 33608.32387 -184010.7447 70785.07178 -34581.30002
 -33608.32387 -184010.7447 70785.07178 -34581.29998
-1.755052267e-13 -167657.9555 109044.9905 -5575.05953
+-1.124859754e-13 -167657.9555 109044.9905 -5575.05953
 65351.59161 -159554.7203 101348.2151 -11312.54113
 93118.40502 -125157.8257 125157.8257 7074.749569
 33822.9268 -139384.3779 139384.3779 17280.58781
@@ -66,7 +66,7 @@
 196824.4702 1.139428721e-12 35498.28088 -55980.07262
 184010.7447 33608.32387 70785.07178 -34581.30005
 184010.7447 -33608.32387 70785.07178 -34581.30004
-167657.9555 1.755052267e-13 109044.9905 -5575.059573
+167657.9555 -1.124859754e-13 109044.9905 -5575.059573
 159554.7203 65351.59161 101348.2151 -11312.54113
 125157.8257 93118.40502 125157.8257 7074.749563
 139384.3779 33822.9268 139384.3779 17280.58778
@@ -117,10 +117,10 @@
 139384.3779 -33822.9268 -139384.3779 17280.58779
 70785.07178 -33608.32387 -184010.7447 46185.62438
 35498.28088 -5.604965797e-12 -196824.4702 53281.72758
-109044.9905 -2.700687402e-13 -167657.9555 36339.18079
+109044.9905 -2.096944977e-13 -167657.9555 36339.18079
 70785.07178 33608.32387 -184010.7447 46185.62438
 184010.7447 -33608.32387 -70785.07178 -34581.30002
-167657.9555 3.016508524e-11 -109044.9905 -5575.059554
+167657.9555 3.006269439e-11 -109044.9905 -5575.059554
 196824.4702 5.258792963e-14 -35498.28088 -55980.07261
 184010.7447 33608.32387 -70785.07178 -34581.30003
 139384.3779 33822.9268 -139384.3779 17280.58778
@@ -130,7 +130,7 @@
 -2.253997394e-13 35498.28088 -196824.4702 53281.72758
 -33608.32387 70785.07178 -184010.7447 46185.62437
 33608.32387 70785.07178 -184010.7447 46185.62438
--2.264523378e-12 109044.9905 -167657.9555 36339.1808
+-2.6687275e-12 109044.9905 -167657.9555 36339.1808
 -65351.59161 101348.2151 -159554.7203 31212.29529
 -93118.40502 125157.8257 -125157.8257 7074.749542
 -33822.9268 139384.3779 -139384.3779 17280.58778
@@ -162,7 +162,7 @@
 -35498.28088 1.269621957e-12 -196824.4702 53281.72758
 -70785.07178 -33608.32387 -184010.7447 46185.62437
 -70785.07178 33608.32387 -184010.7447 46185.62438
--109044.9905 1.048488242e-12 -167657.9555 36339.1808
+-109044.9905 5.609360777e-13 -167657.9555 36339.1808
 -101348.2151 -65351.59161 -159554.7203 31212.29529
 -125157.8257 -93118.40502 -125157.8257 7074.749548
 -139384.3779 -33822.9268 -139384.3779 17280.58776