Merge branch...

Merge branch '140-follow-up-from-wip-resolve-couple-the-transportsimulation-with-richardssimulation' into dg-transport-with-moments

doc/default_files/richards-parameters.xml

@@ -269,7 +269,7 @@ adding an empty line, make text **bold** or ``monospaced``.
     </parameter>
 
     <parameter name="checkTol">
-      <definition> Whenever ``checkJumps`` is activated, it check that flux 
+      <definition> Whenever ``checkJumps`` is activated, it checks that flux 
         reconstruction engine is creating conforming normal fluxes up to ``jumpTol``. 
       </definition>
       <values> float &gt; 0 </values>
@@ -277,6 +277,24 @@ adding an empty line, make text **bold** or ``monospaced``.
     </parameter>
   </category>
 
+  <category name="conservationOfMass">
+    <parameter name="check">
+      <definition> Check that mass is conserved in the domain with respect to
+        boundary fluxes.
+      </definition>
+      <suggestion> none </suggestion>
+      <values> none, warn, error </values>
+    </parameter>
+
+    <parameter name="checkTol">
+      <definition> Whenever ``checkJumps`` is activated, it check that mass is 
+        conserved in the domain with respect to boundary fluxes. ``checkTol`` 
+        is the persentage of mass that can be lost without a warning/error. 
+      </definition>
+      <values> float &gt; 0 </values>
+      <suggestion> 1E-10 </suggestion>
+    </parameter>
+  </category>
 
   <category name="NewtonParameters">
     <parameter name="ForceIteration" hidden="true">

doc/default_files/transport-parameters.xml

@@ -304,7 +304,7 @@ adding an empty line, make text **bold** or ``monospaced``.
     </parameter>
 
     <parameter name="checkTol">
-      <definition> Whenever ``checkJumps`` is activated, it check that flux 
+      <definition> Whenever ``checkJumps`` is activated, it checks that flux 
         reconstruction engine is creating conforming normal fluxes up to ``jumpTol``. 
       </definition>
       <values> float &gt; 0 </values>
@@ -312,6 +312,25 @@ adding an empty line, make text **bold** or ``monospaced``.
     </parameter>
   </category>
 
+  <category name="conservationOfMass">
+    <parameter name="check">
+      <definition> Check that mass is conserved in the domain with respect to
+        boundary fluxes.
+      </definition>
+      <suggestion> none </suggestion>
+      <values> none, warn, error </values>
+    </parameter>
+
+    <parameter name="checkTol">
+      <definition> Whenever ``checkJumps`` is activated, it check that mass is 
+        conserved in the domain with respect to boundary fluxes. ``checkTol`` 
+        is the persentage of mass that can be lost without a warning/error. 
+      </definition>
+      <values> float &gt; 0 </values>
+      <suggestion> 1E-10 </suggestion>
+    </parameter>
+  </category>
+
   <category name="solverParameters">
     <parameter name="reduction">
       <definition> Required relative defect reduction.

dune/dorie/common/flux_reconstruction/raviart_thomas.hh

@@ -283,9 +283,6 @@ public:
   }
 
 private:
-  /// Controls how to deal with jumps occurring in the reconstructed flux
-  enum class CheckLevel { none, warn, error };
-
   std::shared_ptr<spdlog::logger> _log;
   GV _gv;
   FEMU _femu;
@@ -297,25 +294,6 @@ private:
   /// The jump check level
   const CheckLevel _check_lvl;
   const RF _check_tol;
-
-  /// Return a check level for a given string
-  CheckLevel to_check_lvl (const std::string check_lvl_str)
-  {
-    if (to_lower(check_lvl_str) == "none") {
-      return CheckLevel::none;
-    }
-    else if (check_lvl_str == "warn") {
-      return CheckLevel::warn;
-    }
-    else if (check_lvl_str == "error") {
-      return CheckLevel::error;
-    }
-    else {
-      _log->error("Unknown check for jumps in reconstructed flux: {}",
-                  check_lvl_str);
-      DUNE_THROW(IOError, "Unknown check for flux reconstruction jumps");
-    }
-  }
 };
 
 } // namespace Dorie

dune/dorie/common/grid_function_integrator.hh

+#ifndef DUNE_DORIE_GRID_FUNCTION_INTEGRATOR_HH
+#define DUNE_DORIE_GRID_FUNCTION_INTEGRATOR_HH
+
+namespace Dune {
+namespace Dorie {
+
+/*-------------------------------------------------------------------------*//**
+ * @brief      Integrate grid function
+ * 
+ * @author     Santiago Ospina De Los Ríos
+ * @date       2019
+ * 
+ * @param[in]  gf         The grid function
+ * @param[in]  int_order  The integration order
+ *
+ * @tparam     GF         The grid function type
+ *
+ * @return     The integral over the whole domain
+ */
+template<class GF>
+typename GF::Traits::RangeType
+gf_integrator(const GF& gf, unsigned int int_order)
+{
+  using Traits = typename GF::Traits;
+  using Range = typename Traits::RangeType;
+
+  Range integral(0.);
+
+  const auto& gv = gf.getGridView();
+
+  // Loop over the grid
+  for (const auto& cell : Dune::elements(gv,Dune::Partitions::interior))
+  {
+    // Get geometry
+    auto geo = cell.geometry();
+
+    // Integrate on the cell
+    for (const auto& point : quadratureRule(geo, int_order))
+    {
+      // Evaluate function
+      Range y;
+      const auto& x = point.position();
+      gf.evaluate(cell,x,y);
+
+      auto factor = point.weight() * geo.integrationElement(x);
+
+      // Accumulate
+      integral += y*factor;
+    }
+  }
+
+  // Accumulate over processors
+  return gv.comm().sum(integral);
+}
+
+/*-------------------------------------------------------------------------*//**
+ * @brief      Integrate grid function on the boundary of the domain
+ * @details    If the grid function is a vector valued function, it will be
+ *             integrated in the normal direction of the boundary.
+ * 
+ * @author     Santiago Ospina De Los Ríos
+ * @date       2019
+ * 
+ * @param[in]  gf         The grid function
+ * @param[in]  int_order  The integration order
+ *
+ * @tparam     GF         The grid function type
+ *
+ * @return     The integral over the whole domain boundary
+ */
+template<class GF>
+Dune::FieldVector<typename GF::Traits::RangeFieldType,1>
+boundary_gf_integrator(const GF& gf, unsigned int int_order)
+{
+  using Traits = typename GF::Traits;
+  using Range = typename Traits::RangeType;
+  constexpr int dim_range = Traits::dimRange;
+  constexpr int dim_domain = Traits::dimDomain;
+
+  static_assert(dim_range == 1 || dim_range == dim_domain);
+
+  Dune::FieldVector<typename GF::Traits::RangeFieldType,1> integral(0.);
+
+  const auto& gv = gf.getGridView();
+
+  // Loop over the grid
+  for (const auto& cell : Dune::elements(gv,Dune::Partitions::interior))
+  {
+    // Loop over the intersections
+    for (const auto& face : Dune::intersections(gv,cell)) 
+    {
+      // Only integrate on boundary
+      if (face.boundary() and not face.neighbor())
+      {
+        // Get geometry
+        auto geo_f = face.geometry();
+        auto geo_in_i = face.geometryInInside();
+
+        // Integrate on the face
+        for (const auto& point : quadratureRule(geo_f, int_order))
+        {
+          // Relative positions
+          const auto& position_f = point.position();
+          const auto position_i = geo_in_i.global(position_f);
+
+          // Face normal vector
+          const auto normal = face.unitOuterNormal(position_f);
+
+          // Evaluate function
+          Range y;
+          gf.evaluate(cell,position_i,y);
+
+          const auto factor = point.weight() * geo_f.integrationElement(position_f);
+          
+          // Accumulate
+          if constexpr (dim_range == dim_domain)
+            integral += y * normal * factor;
+          else 
+            integral += y * factor;
+        }
+      }
+    }
+  }
+
+  // Accumulate over processors
+  return gv.comm().sum(integral);
+}
+
+}
+}
+
+#endif

dune/dorie/common/statistical_moments.hh

+#ifndef DUNE_DORIE_STATISTICAL_MOMENTS_HH
+#define DUNE_DORIE_STATISTICAL_MOMENTS_HH
+
+namespace Dune {
+namespace Dorie {
+
+template<class GF>
+std::vector<typename GF::Traits::DomainType> 
+statistical_moments(
+  const GF& gf, 
+  unsigned int int_order, 
+  unsigned int max_order = 2)
+{
+  using Traits = typename GF::Traits;
+  using Domain = typename Traits::DomainType;
+  constexpr unsigned int dim = Traits::dimDomain;
+
+  assert(max_order>0);
+
+  std::vector<Domain> moment(max_order);
+
+  const auto& gv = gf.getGridView();
+
+  // Loop over the grid
+  for (const auto& cell : Dune::elements(gv,Dune::Partitions::interior))
+  {
+    auto geo = cell.geometry();
+    for (const auto& point : quadratureRule(geo, int_order))
+    {
+      const auto& x = point.position();
+      auto x_gl = geo.global(x);
+
+      typename Traits::RangeType y;
+      gf.evaluate(cell,x,y);
+
+      auto factor = point.weight() * geo.integrationElement(x);
+
+      Domain weight(1.);
+      for (unsigned int i = 0; i < max_order; ++i) {
+        for (unsigned int j = 0; j < dim; ++j)
+        {
+          weight[j] *= x_gl[j];
+          moment[i][j] += weight[j]*y*factor;
+        }
+      }
+      
+    }
+  }
+  
+  // Accumulate over processors
+  for (unsigned int i = 0; i < max_order; ++i)
+    moment[i] = gv.comm().sum(moment[i]);
+  return moment;
+}
+
+}
+}
+
+#endif

dune/dorie/common/util.hh

@@ -31,6 +31,10 @@ enum class OutputPolicy {None,EndOfRichardsStep,EndOfTransportStep};
  */
 enum class AdaptivityPolicy {None,WaterFlux,SoluteFlux};
 
+/// Controls how to deal with run-time checks
+enum class CheckLevel { none, warn, error };
+
+
 /// Return the estimation of entries per matrix row for the spatial GridOperator
 /** This supposedly decreases matrix assembly time.
  *  The values specify the *blocks* per row. DG assembles one block for the

dune/dorie/common/utility.hh

@@ -14,6 +14,7 @@
 #include <dune/common/fvector.hh>
 
 #include <dune/grid/common/gridview.hh>
+#include <dune/dorie/common/util.hh>
 
 namespace Dune {
 namespace Dorie {
@@ -116,6 +117,23 @@ inline bool is_none (const std::string& in)
     return false;
 }
 
+/// Return a check level for a given string
+inline CheckLevel to_check_lvl (const std::string check_lvl_str)
+{
+  if (to_lower(check_lvl_str) == "none") {
+    return CheckLevel::none;
+  }
+  else if (check_lvl_str == "warn") {
+    return CheckLevel::warn;
+  }
+  else if (check_lvl_str == "error") {
+    return CheckLevel::error;
+  }
+  else {
+    DUNE_THROW(IOError, "Unknown check level");
+  }
+}
+
 /*-------------------------------------------------------------------------*//**
  * @brief      Transform a Yaml node with sequences to map
  * @details    This function check all the nodes within 'in' and transform every

dune/dorie/model/richards/richards.cc

+#include <dune/pdelab/function/scalarscaled.hh>
 #include <dune/dorie/model/richards/richards.hh>
+#include <dune/dorie/common/grid_function_integrator.hh>
+#include <dune/dorie/common/statistical_moments.hh>
+
+#include <limits>
+#include <cmath>
 
 namespace Dune{
 namespace Dorie{
@@ -94,6 +100,16 @@ RichardsSimulation<Traits>::RichardsSimulation (
 							inifile.get<std::string>("output.fileName"),
 							inifile.get<std::string>("output.outputPath"),
 							"./");
+
+    // Configure statistics file
+    std::ofstream stats_file;
+    stats_file.open(inifile.get<std::string>("output.fileName") + ".dat", std::ofstream::trunc);
+    stats_file << "# Statistics of the Richards Equation generated by DORiE.\n";
+    stats_file << "#\n";
+    stats_file << "# run_time\t time\t mass\t boundary_flux\t";
+    stats_file << "mean_cw_x\t mean_cw_y\t mean_cw_z\t";
+    stats_file << "variance_cw_x\t variance_cw_y\t variance_cw_z\n";
+    stats_file.close();
 	}
 
 	// --- Adaptivity Setup --- //
@@ -229,6 +245,8 @@ void RichardsSimulation<Traits>::step()
 	// invalidate cache for water flux reconstruction
 	cache_water_flux_gf_rt.reset();
 
+	time_steps++;
+
 	if (this->output_policy() == OutputPolicy::EndOfRichardsStep)
 		write_data();
 }
@@ -266,15 +284,102 @@ void RichardsSimulation<Traits>::post_adapt_grid()
 template<typename Traits>
 void RichardsSimulation<Traits>::write_data () const
 {
-	if (vtkwriter)
+	const auto time = controller->getTime();
+	
+	// create adapters
+	auto head = get_matric_head();
+	auto wflux = std::make_shared<const GFWaterFlux>(gfs, u, fparam);
+	auto cond = std::make_shared<const GFConductivity>(gv, fparam);
+	auto wc =  std::make_shared<const GFWaterContent>(head, gv, fparam);
+	auto sat = std::make_shared<const GFSaturation>(head, gv, fparam);
+
+	// Compute mass of WC
+	RF mass = gf_integrator(*wc,2*order+2);
+
+	// Compute statistical moments of WC
+	using PDF = Dune::PDELab::ScalarScaledGridFunctionAdapter<const GFWaterContent>;
+	PDF pdf(1./mass,*wc);
+	auto moments = statistical_moments(pdf,order*4,2);
+	auto moments_squared = moments[0];
+	for (int i = 0; i < dim; ++i) moments_squared[i]*=moments_squared[i];
+	const auto variance = moments[1] - moments_squared;
+	Dune::FieldVector<double,3> moments0_3d(0.), variance_3d(0.);
+	for (int i = 0; i < dim; ++i)
 	{
-		// create adapters
-		auto head = get_matric_head();
-		auto wflux = std::make_shared<const GFWaterFlux>(gfs, u, fparam);
-		auto cond = std::make_shared<const GFConductivity>(gv, fparam);
-		auto wc =  std::make_shared<const GFWaterContent>(head, gv, fparam);
-		auto sat = std::make_shared<const GFSaturation>(head, gv, fparam);
+		moments0_3d[i]  = moments[0][i];
+		variance_3d[i]  = variance[i];
+	}
+
+	RF boundary_flux = std::numeric_limits<RF>::quiet_NaN();;
+	if constexpr (enable_rt_engine) {
+		if (enable_fluxrc) {
+			auto wfluxr = get_water_flux_reconstructed();
+		    boundary_flux = boundary_gf_integrator(*wfluxr,2*flux_order+2);
+		}
+	}
+
+    if (time_steps==0) cache_time          = time;
+    if (time_steps==0) cache_mass          = mass;
+    if (time_steps==0) cache_boundary_flux = boundary_flux;
 
+	// hard check on mass
+	if (std::isnan(mass)) {
+		this->_log->error("Mass of water is NaN: PDE solver has diverged.");
+		DUNE_THROW(MathError,"PDE solver has diverged!");
+	}
+    // check relative change of mass
+    if (Dune::FloatCmp::ne(mass,0.))
+    {
+		// change in time
+		auto dt = time - cache_time;
+
+		// change of mass by PDE solution
+		auto mass_dt = (mass - cache_mass);
+
+		// change of mass by boundary conditions
+		mass_dt += 0.5*dt*(boundary_flux + cache_boundary_flux);
+
+		// relative mass change
+		auto rel_mass_dt = mass_dt/mass;
+		auto mass_tol = inifile.get<double>("conservationOfMass.checkTol");
+		auto check_mass_str = inifile.get<std::string>("conservationOfMass.check");
+		auto check_mass_lvl = Dorie::to_check_lvl(check_mass_str);
+
+		const bool passed = rel_mass_dt < mass_tol;
+
+		if (passed) {
+		  this->_log->trace("  Conservation of mass does not hold for"
+		                    " {}% of mass",
+		                    rel_mass_dt);
+		}
+		else if (not passed && check_mass_lvl == CheckLevel::warn) {
+		  this->_log->warn("Conservation of mass above tolerance:"
+		                    " {}% of mass is not conserved.",
+		                    rel_mass_dt);
+		}
+		else if (not passed && check_mass_lvl == CheckLevel::error) {
+		  this->_log->error("Conservation of mass above tolerance:"
+		                    " {}% of mass is not conserved.",
+		                    rel_mass_dt);
+		  DUNE_THROW(MathError, "Conservation of mass does not obey tolerance");
+		}
+    }
+
+    // update cache
+    cache_time          = time;
+    cache_mass          = mass;
+    cache_boundary_flux = boundary_flux;
+
+	// Configure statistics file
+	std::ofstream stats_file;
+	stats_file.open(inifile.get<std::string>("output.fileName") + ".dat", std::ofstream::app);
+	stats_file << timer.elapsed() << "\t" << time << "\t" << mass << "\t";
+	stats_file << boundary_flux << "\t" << moments0_3d << "\t" << variance_3d;
+	stats_file << std::endl;
+	stats_file.close();
+
+	if (vtkwriter)
+	{
 		if (inifile.get<bool>("output.vertexData")) {
 			vtkwriter->addVertexData(head,"head");
 			vtkwriter->addVertexData(wflux,"flux");
@@ -302,7 +407,6 @@ void RichardsSimulation<Traits>::write_data () const
 		}
 
 		try{
-			const auto time = controller->getTime();
 			this->_log->trace("Writing solution at time {:.2e}", time);
 			vtkwriter->write(time, output_type);
 		}

dune/dorie/model/richards/richards.hh

@@ -15,6 +15,7 @@
 
 #include <dune/dorie/common/simulation_base.hh>
 #include <dune/dorie/common/util.hh>
+#include <dune/dorie/common/utility.hh>
 #include <dune/dorie/common/interpolator.hh>
 #include <dune/dorie/common/time_controller.hh>
 #include <dune/dorie/common/grid_creator.hh>
@@ -190,6 +191,8 @@ private:
 	static constexpr int order = Traits::order;
 	static constexpr int flux_order = Traits::flux_order;
 	static constexpr bool enable_rt_engine = Traits::enable_rt_engine;
+	
+	using TimeField = typename Traits::TimeField;
 	using RF = typename Traits::RF;
 	using Grid = typename Traits::Grid;
 	using GV = typename Traits::GV;
@@ -316,12 +319,18 @@ protected:
 	std::unique_ptr<Writer> vtkwriter;
 	std::unique_ptr<AdaptivityHandler> adaptivity;
 
-	mutable std::shared_ptr<GFFluxReconstruction> cache_water_flux_gf_rt;
+	Dune::Timer timer;
 
 	RF time_before;
 	RF dt_before;
-
 	const bool enable_fluxrc;
+  	unsigned int time_steps;
+  	
+
+	mutable std::shared_ptr<GFFluxReconstruction> cache_water_flux_gf_rt;  	
+  	mutable RF cache_mass;
+  	mutable RF cache_boundary_flux;
+  	mutable TimeField cache_time;
 public:
 
 	/*------------------------------------------------------------------------*//**

dune/dorie/model/transport/transport.cc

+#include <dune/pdelab/function/scalarscaled.hh>
 #include <dune/dorie/model/transport/transport.hh>
+#include <dune/dorie/common/grid_function_integrator.hh>
+#include <dune/dorie/common/statistical_moments.hh>
+
+
+#include <limits>
+#include <cmath>
 
 namespace Dune{
 namespace Dorie{
@@ -85,6 +92,16 @@ TransportSimulation<Traits>::TransportSimulation(
               inifile.get<std::string>("output.fileName"),
               inifile.get<std::string>("output.outputPath"),
               "./");
+
+    // Configure statistics file
+    std::ofstream stats_file;
+    stats_file.open(inifile.get<std::string>("output.fileName") + ".dat", std::ofstream::trunc);
+    stats_file << "# Statistics of the solute transport equation generated by DORiE.\n";
+    stats_file << "#\n";
+    stats_file << "# run_time\t time\t mass\t boundary_flux\t";
+    stats_file << "mean_cw_x\t mean_cw_y\t mean_cw_z\t";
+    stats_file << "variance_cw_x\t variance_cw_y\t variance_cw_z\n";
+    stats_file.close();
   }
 
   // --- Operator Setup ---
@@ -277,11 +294,100 @@ void TransportSimulation<Traits>::step()
 template<typename Traits>
 void TransportSimulation<Traits>::write_data () const
 {
+  const auto time = controller->getTime();
+
+  auto solute_gf = std::make_shared<GFSolute>(gfs,u);
+  auto solute_t_gf = std::make_shared<GFTotalSolute>(solute_gf, water_content_gf);
+
+  // Compute statistical moments of Cw
+  using PDF = Dune::PDELab::ScalarScaledGridFunctionAdapter<GFSolute>;
+  PDF pdf(1./gf_integrator(*solute_gf,2*order+2),*solute_gf);
+  auto moments = statistical_moments(pdf,order*4,2);
+  auto moments_squared = moments[0];
+  for (int i = 0; i < dim; ++i) moments_squared[i]*=moments_squared[i];
+  const auto variance = moments[1] - moments_squared;
+  Dune::FieldVector<double,3> moments0_3d(0.), variance_3d(0.);
+  for (int i = 0; i < dim; ++i)
+  {
+    moments0_3d[i]  = moments[0][i];
+    variance_3d[i]  = variance[i];
+  }
+
+  RF mass = gf_integrator(*solute_t_gf,2*order+2);
+  RF boundary_flux = std::numeric_limits<RF>::quiet_NaN();
+  if constexpr (enable_rt_engine) {
+    if (enable_fluxrc) {
+      auto cwfluxr = get_solute_flux_reconstructed();
+      boundary_flux = boundary_gf_integrator(*cwfluxr,2*flux_order+2);
+    }
+  }
+
+  if (time_steps==0) cache_time          = time;
+  if (time_steps==0) cache_mass          = mass;
+  if (time_steps==0) cache_boundary_flux = boundary_flux;
+
+  // hard check on mass
+  if (std::isnan(mass)) {
+    this->_log->error("Mass of solute is NaN: PDE solver has diverged.");
+    DUNE_THROW(MathError,"PDE solver has diverged!");
+  }
+
+  // check relative change of mass
+  if (Dune::FloatCmp::ne(mass,0.))
+  {
+    // change in time
+    auto dt = time - cache_time;
+
+    // change of mass by PDE solution
+    auto mass_dt = (mass - cache_mass);
+
+    // change of mass by boundary conditions
+    mass_dt += 0.5*dt*(boundary_flux + cache_boundary_flux);
+
+    // relative mass change
+    auto rel_mass_dt = mass_dt/mass;
+    auto mass_tol = inifile.get<double>("conservationOfMass.checkTol");
+    auto check_mass_str = inifile.get<std::string>("conservationOfMass.check");
+    auto check_mass_lvl = Dorie::to_check_lvl(check_mass_str);
+
+    const bool passed = rel_mass_dt < mass_tol;
+
+    if (passed) {
+      this->_log->trace("  Conservation of mass does not hold for"
+                        " {}% of mass",
+                        rel_mass_dt);
+    }
+    else if (not passed && check_mass_lvl == CheckLevel::warn) {
+      this->_log->warn("Conservation of mass above tolerance:"