Support plotting fluxes through KnoFuPlotter

Distinguish plotting of states and fluxes. A separate call to plot
fluxes now reconstructs fluxes and plots them as vertex data into the
specified path.

dune/dorie/model/richards/knofu/plotting.hh

@@ -28,15 +28,18 @@ private:
     using U = typename Base::U;
     /// Type of the discrete grid function based on the solution vector
     using GFSolution = typename Base::GFSolution;
+    /// Type of the flux grid function
+    using GFFluxReconstruction = typename Traits::GFFluxReconstruction;
     /// Type of the VTK function adapter based on the solution grid function
-    using VTKFunction = Dune::PDELab::VTKGridFunctionAdapter<GFSolution>;
+    template<class GF>
+    using VTKFunction = Dune::PDELab::VTKGridFunctionAdapter<GF>;
 
     /// VTK writer
     Dune::SubsamplingVTKWriter<typename Traits::GV> vtkwriter;
-    /// State storage
-    std::vector<std::shared_ptr<U>> member_states;
-    /// Name storage associated with respective state
-    std::vector<std::string> member_names;
+    /// State storage (name, solution, time)
+    std::vector<std::tuple<std::string,
+                           std::shared_ptr<U>,
+                           double>> member_storage;
     /// Output file encoding
     const Dune::VTK::OutputType output_type;
 
@@ -66,46 +69,114 @@ public:
      *  supplied.
      *  \param state State vector (water content)
      *  \param name Name of the state (i.e. member number)
+     *  \param time The simulation time associated with this state
+     *              (required for flux reconstruction)
      */
     void add_solution_to_output (
         const BackendVector& state,
-        const std::string& name)
+        const std::string& name,
+        const double time)
     {
-        auto solution = std::make_shared<U>(*this->gfs, 0.0);
-        Dune::PDELab::Backend::native(*solution) = state;
+        auto state_copy = std::make_shared<U>(*this->gfs, 0.0);
+        Dune::PDELab::Backend::native(*state_copy) = state;
 
-        member_states.push_back(solution);
-        member_names.push_back(name);
+        member_storage.emplace_back(name, state_copy, time);
     }
 
-    /// Write all stored states into a .vtu file. Then clear all data
+    /// Release the data stored in this object
+    void clear ()
+    {
+        // release data
+        member_storage.clear();
+    }
+
+    /// Write all stored states into a .vtu file.
     /** Create discrete grid functions and VTK functions from the saved
      *  states.
      *  \param path Path to write into
      *  \param name Filename to write
      */
-    void write_output (const std::string& path, const std::string& name)
+    void write_state (const std::string& path, const std::string& name)
+    {
+        const auto vertex_data
+            = this->inifile.template get<bool>("output.vertexData");
+
+        for(auto&& member : member_storage)
+        {
+            // fetch data from storage
+            const auto& state = std::get<std::shared_ptr<U>>(member);
+            const auto& name = std::get<std::string>(member);
+
+            // stack VTKFunctions into writer
+            auto solution = std::make_shared<GFSolution>(this->gfs,
+                                                         state);
+            using VTKF = VTKFunction<GFSolution>;
+            add_vtk_data(std::make_shared<VTKF>(solution, name), vertex_data);
+        }
+
+        write_vtk_output(path, name);
+    }
+
+    /// Transform all stored states into fluxes and write these into a file
+    void write_flux (const std::string& path, const std::string& name)
     {
-        // stack VTKFunctions into writer
-        for(size_t i = 0; i < member_states.size(); ++i) {
-            auto dgf = std::make_shared<GFSolution>(*this->gfs,
-                                                    *member_states.at(i));
-            auto dgf_vtk = std::make_shared<VTKFunction>(dgf,
-                                                         member_names.at(i));
-
-            if (this->inifile.template get<bool>("output.vertexData"))
-                vtkwriter.addVertexData(dgf_vtk);
-            else
-                vtkwriter.addCellData(dgf_vtk);
+        for(auto&& member : member_storage)
+        {
+            // fetch data from storage
+            const auto& state = std::get<std::shared_ptr<U>>(member);
+            const auto& name = std::get<std::string>(member);
+            const auto& time = std::get<double>(member);
+
+            // copy the state
+            const auto state_head = std::make_shared<U>(*state);
+
+            // Transform water content into matric head
+            using GFWaterContentToHead = typename Base::GFWaterContentToHead;
+            if (this->_solution_type == FilterSolutionType::water_content)
+            {
+                const auto wc_solution = std::make_shared<GFSolution>(this->gfs,
+                                                                      state);
+                GFWaterContentToHead head_gf(this->inifile,
+                                             this->gv,
+                                             wc_solution,
+                                             this->fparam);
+                Dune::PDELab::interpolate(head_gf, *this->gfs, *state_head);
+            }
+
+            // reconstruct the flux
+            typename Traits::State state_pack{this->gfs, state_head, time};
+            auto dgf = this->get_water_flux_reconstructed(state_pack);
+
+            // stack VTKFunctions into writer
+            using VTKF = VTKFunction<GFFluxReconstruction>;
+            add_vtk_data(std::make_shared<VTKF>(dgf, name), true);
         }
 
-        // write output
+        write_vtk_output(path, name);
+    }
+
+private:
+
+    /// Write the data stored in the VTK writer
+    void write_vtk_output (const std::string& path,
+                           const std::string& name,
+                           const bool clear_writer=true)
+    {
         vtkwriter.pwrite(name, path, "", output_type);
 
-        // release data
-        vtkwriter.clear();
-        member_states.clear();
-        member_names.clear();
+        if (clear_writer)
+            vtkwriter.clear();
+    }
+
+    /// Add a VTK grid function to the
+    template<typename VTKGridFunction>
+    void add_vtk_data (const std::shared_ptr<VTKGridFunction> vtk_gf,
+                       const bool vertex_data)
+    {
+        if (vertex_data)
+            vtkwriter.addVertexData(vtk_gf);
+        else
+            vtkwriter.addCellData(vtk_gf);
     }
 
 };

[Santiago Ospina De Los Ríos]

@lriedel Remember that the flux reconstruction engine is not available for every polynomial order. Use the enable_rt_engine bool in Traits to either ensure that this will always work (e.g. not printing when unavailable), or by restricting KnoFu interface for valid combinations of the engine.

[Lukas Riedel]

Thanks for the heads-up! We currently do the latter: The interface is restricted to low-order cases, namely orders 0 and 1 for now.

[Santiago Ospina De Los Ríos]

Yeah, I know, I am just suggesting a static_assert so that higher-order instantiations (who knows) have an appropriate error message and not a messy and long template error message.