Merge pull request geodynamics#2851 from gassmoeller/add_option_for_h…

…igher_order_output

Add option for higher order output

doc/modules/changes/20190313_gassmoeller_2

@@ -0,0 +1,8 @@
+New: There is now an option to output visualization data as higher order
+polynomials. This is an improvement in accuracy and requires less disk space
+than the 'Interpolate output' option that was available before (the new option
+requires the old option to be set, which it is by default). However the new
+output can only be read by ParaView version 5.5 and newer and is therefore
+disabled by default.
+<br>
+(Rene Gassmoeller, 2019/03/13)

include/aspect/postprocess/visualization.h

@@ -462,6 +462,19 @@ namespace aspect
          */
         bool filter_output;
 
+        /**
+         * deal.II offers the possibility to write vtu files with higher order
+         * representations of the output data. This means each cell will correctly
+         * show the higher order representation of the output data instead of the
+         * linear interpolation between vertices that ParaView and Visit usually show.
+         * Note that activating this option is safe and recommended, but requires that
+         * (i) ``Output format'' is set to ``vtu'', (ii) ``Interpolate output'' is
+         * set to true, (iii) you use a sufficiently new version of Paraview
+         * or Visit to read the files (Paraview version 5.5 or newer, and Visit version
+         * to be determined), and (iv) you use deal.II version 9.1.0 or newer.
+         */
+        bool write_higher_order_output;
+
         /**
          * For free surface computations Aspect uses an Arbitrary-Lagrangian-
          * Eulerian formulation to handle deforming the domain, so the mesh

source/postprocess/visualization.cc

@@ -573,6 +573,10 @@ namespace aspect
           vtk_flags.cycle = this->get_timestep_number();
           vtk_flags.time = time_in_years_or_seconds;
 
+#if DEAL_II_VERSION_GTE(9,1,0)
+          vtk_flags.write_higher_order_cells = write_higher_order_output;
+#endif
+
           data_out.set_flags (vtk_flags);
 
           // Write as many files as processes. For this case we support writing in a
@@ -838,6 +842,30 @@ namespace aspect
                              "and a factor of 8 in 3d, when using quadratic elements for the velocity, "
                              "and correspondingly more for even higher order elements.");
 
+          prm.declare_entry ("Write higher order output", "false",
+                             Patterns::Bool(),
+                             "deal.II offers the possibility to write vtu files with higher order "
+                             "representations of the output data. This means each cell will correctly "
+                             "show the higher order representation of the output data instead of the "
+                             "linear interpolation between vertices that ParaView and Visit usually show. "
+                             "Note that activating this option is safe and recommended, but requires that "
+                             "(i) ``Output format'' is set to ``vtu'', (ii) ``Interpolate output'' is "
+                             "set to true, (iii) you use a sufficiently new version of Paraview "
+                             "or Visit to read the files (Paraview version 5.5 or newer, and Visit version "
+                             "to be determined), and (iv) you use deal.II version 9.1.0 or newer. "
+                             "\n"
+                             "The effect of using this option can be seen in the following "
+                             "picture:"
+                             "\n\n"
+                             "\\begin{center}"
+                             "  \\includegraphics[width=0.5\\textwidth]{viz/parameters/higher-order-output}"
+                             "\\end{center}"
+                             "The top figure shows the plain output without interpolation or higher "
+                             "order output. The middle figure shows output that was interpolated as "
+                             "discussed for the ``Interpolate output'' option. The bottom panel "
+                             "shows higher order output that achieves better accuracy than the "
+                             "interpolated output at a lower memory cost.");
+
           prm.declare_entry ("Filter output", "false",
                              Patterns::Bool(),
                              "deal.II offers the possibility to filter duplicate vertices for HDF5 "
@@ -949,6 +977,25 @@ namespace aspect
 
           interpolate_output = prm.get_bool("Interpolate output");
           filter_output = prm.get_bool("Filter output");
+          write_higher_order_output = prm.get_bool("Write higher order output");
+
+#if DEAL_II_VERSION_GTE(9,1,0)
+#else
+          AssertThrow(write_higher_order_output == false, ExcMessage("The 'Write higher order output' functionality is only "
+                                                                     "available for deal.II version 9.1.0 or newer. Please update "
+                                                                     "your deal.II version if you need this option."));
+#endif
+
+          if (write_higher_order_output == true)
+            {
+              AssertThrow(output_format == "vtu",
+                          ExcMessage("The option 'Postprocess/Visualization/Write higher order output' requires the "
+                                     "data output format to be set to 'vtu'"));
+              AssertThrow(interpolate_output == true,
+                          ExcMessage("The input parameter 'Postprocess/Visualization/Write higher order output' "
+                                     "requires the input parameter 'Interpolate output' to be set to 'true'"));
+            }
+
           output_mesh_velocity = prm.get_bool("Output mesh velocity");
 
           // now also see which derived quantities we are to compute

tests/visualization_higher_order_output.prm.deal.II.9.1

@@ -0,0 +1,17 @@
+set Dimension = 2
+
+include $ASPECT_SOURCE_DIR/cookbooks/van-keken-smooth.prm
+
+set End time = 0.0
+
+subsection Mesh refinement
+  set Initial adaptive refinement              = 0
+  set Initial global refinement                = 3
+end
+
+subsection Postprocess
+  subsection Visualization
+    set Interpolate output = true
+    set Write higher order output = true
+  end
+end

tests/visualization_higher_order_output/screen-output

@@ -0,0 +1,22 @@
+
+Number of active cells: 64 (on 4 levels)
+Number of degrees of freedom: 1,237 (578+81+289+289)
+
+*** Timestep 0:  t=0 seconds
+   Skipping temperature solve because RHS is zero.
+   Solving C_1 system ... 0 iterations.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system... 18+0 iterations.
+
+   Postprocessing:
+     Writing graphical output:  output-visualization_higher_order_output/solution/solution-00000
+     RMS, max velocity:         0.000249 m/s, 0.000589 m/s
+     Compositions min/max/mass: 0/1/0.1874
+
+Termination requested by criterion: end time
+
+
++----------------------------------------------+------------+------------+
++----------------------------------+-----------+------------+------------+
++----------------------------------+-----------+------------+------------+
+

tests/visualization_higher_order_output/statistics

@@ -0,0 +1,19 @@
+# 1: Time step number
+# 2: Time (seconds)
+# 3: Time step size (seconds)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Iterations for temperature solver
+# 9: Iterations for composition solver 1
+# 10: Iterations for Stokes solver
+# 11: Velocity iterations in Stokes preconditioner
+# 12: Schur complement iterations in Stokes preconditioner
+# 13: Visualization file name
+# 14: RMS velocity (m/s)
+# 15: Max. velocity (m/s)
+# 16: Minimal value for composition C_1
+# 17: Maximal value for composition C_1
+# 18: Global mass for composition C_1
+0 0.000000000000e+00 0.000000000000e+00 64 659 289 289 0 0 17 19 18 output-visualization_higher_order_output/solution/solution-00000 2.48863895e-04 5.88814545e-04 0.00000000e+00 1.00000000e+00 1.87415956e-01