add alpha_boundary function to error indicator

FlowBoundary and time now need to be passed to the error indicator.

dune/dorie/interface/adaptivity.hh

@@ -5,12 +5,13 @@ namespace Dune{
 namespace Dorie{
 
 /// Adaptivity base class. Does nothing.
-template<typename Traits, typename GFS, typename Param, typename U>
+template<typename Traits, typename GFS, typename Param, typename Boundary, typename U>
 class AdaptivityHandlerBase
 {
 private:
 	using Grid = typename Traits::Grid;
 	using GV = typename Traits::GV;
+	using RF = typename Traits::RF;
 
 public:
 	AdaptivityHandlerBase () = default;
@@ -19,12 +20,16 @@ public:
 	/// Do nothing.
 	/** \return Boolean if operators have to be updated
 	 */
-	virtual bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, Param& param, U& uold, U& unew) { return false; }
+	virtual bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& uold, U& unew) { return false; }
 };
 
 /// Specialized SimulationBase class providing functions and members for adaptive grid refinement
-template<typename Traits, typename GFS, typename Param, typename U>
-class AdaptivityHandler : public AdaptivityHandlerBase<Traits,GFS,Param,U>
+template<typename Traits,
+	typename GFS,
+	typename Param,
+	typename Boundary,
+	typename U>
+class AdaptivityHandler : public AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U>
 {
 private:
 	static constexpr int order = Traits::fem_order;
@@ -37,7 +42,7 @@ private:
 	/// Adaptivity GFS
 	using AGFS = typename AGFSHelper::Type;
 	/// Error estimator local operator
-	using ESTLOP = Dune::Dorie::FluxErrorEstimator<Traits,Param,typename AGFSHelper::FEM>;
+	using ESTLOP = Dune::Dorie::FluxErrorEstimator<Traits,Param,Boundary>;
 	/// Empty constraints container
 	using NoTrafo = Dune::PDELab::EmptyTransformation;
 	/// Solution vector type
@@ -66,7 +71,7 @@ public:
 	 *  \param grid Grid to adapt (reference is not saved)
 	 */
 	AdaptivityHandler (Dune::ParameterTree& _inifile, Grid& grid)
-		: AdaptivityHandlerBase<Traits,GFS,Param,U>(),
+		: AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U>(),
 		inifile(_inifile),
 		maxLevel(inifile.get<int>("adaptivity.maxLevel")),
 		minLevel(inifile.get<int>("adaptivity.minLevel")),
@@ -95,7 +100,15 @@ public:
 	 *  \param unew Solution vector
 	 *  \return Boolean if grid operators have to be set up again (true if grid changed)
 	 */
-	bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, Param& param, U& uold, U& unew) override
+	bool adapt_grid (
+		Grid& grid,
+		GV& gv,
+		GFS& gfs,
+		const Param& param,
+		const Boundary& boundary,
+		const RF time,
+		U& uold,
+		U& unew) override
 	{
 		Dune::Timer timer2;
 		float t2;
@@ -118,7 +131,8 @@ public:
 
 			// set up helper GFS
 			AGFS p0gfs = AGFSHelper::create(gv);
-			ESTLOP estlop(inifile,param);
+			ESTLOP estlop(inifile, param, boundary);
+			estlop.setTime(time);
 			MBE mbe(0);
 			ESTGO estgo(gfs,p0gfs,estlop,mbe);
 			U0 eta(p0gfs,0.0);
@@ -205,14 +219,14 @@ public:
  *  The constexpr boolean adapt_grid in Dune::Dorie::BaseTraits will
  *  define which create() function is enabled at compile time.
  */
-template<typename Traits, typename GFS, typename Param, typename U>
+template<typename Traits, typename GFS, typename Param, typename Boundary, typename U>
 struct AdaptivityHandlerFactory
 {
 private:
 	static constexpr bool enabled = Traits::adapt_grid;
 	using Grid = typename Traits::Grid;
 
-	using AHB = AdaptivityHandlerBase<Traits,GFS,Param,U>;
+	using AHB = AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U>;
 
 	Dune::ParameterTree& inifile;
 	Grid& grid;
@@ -246,7 +260,7 @@ public:
 									std::is_same<Grid,Dune::UGGrid<3>>::value,
 									"Adaptivity is only implemented for UGGrid!");
 		std::unique_ptr<AHB> p;
-		p = std::make_unique<AdaptivityHandler<Traits,GFS,Param,U>>(inifile,grid);
+		p = std::make_unique<AdaptivityHandler<Traits,GFS,Param,Boundary,U>>(inifile,grid);
 		return p;
 	}
 };

dune/dorie/interface/simulation.cc

@@ -165,12 +165,13 @@ bool Simulation<Traits>::compute_time_step ()
 template<typename Traits>
 void Simulation<Traits>::run ()
 {
-	output->write_vtk_output(controller->getTime());
+	const auto t_start = controller->getTime();
+	output->write_vtk_output(t_start);
 	while(controller->doStep()) {
 		if(!compute_time_step()){
 			continue;
 		}
-		if(adaptivity->adapt_grid(*grid, gv, *gfs, *param, *uold, *unew)){ // reset operators if grid changes
+		if(adaptivity->adapt_grid(*grid, gv, *gfs, *param, *fboundary, t_start, *uold, *unew)){ // reset operators if grid changes
 			operator_setup();
 		}
 		output->write_vtk_output(controller->getTime());

dune/dorie/interface/simulation.hh

@@ -83,9 +83,9 @@ protected:
 	/// Factory for creating the output writer
 	using OutputWriterFactory = Dune::Dorie::OutputWriterFactory<Traits,GFS,Parameters,U>;
 	/// Grid Adaptivity Handler base class
-	using AdaptivityHandler = Dune::Dorie::AdaptivityHandlerBase<Traits,GFS,Parameters,U>;
+	using AdaptivityHandler = Dune::Dorie::AdaptivityHandlerBase<Traits,GFS,Parameters,FlowBoundary,U>;
 	/// Factory for creating the AdaptivityHandler
-	using AdaptivityHandlerFactory = Dune::Dorie::AdaptivityHandlerFactory<Traits,GFS,Parameters,U>;
+	using AdaptivityHandlerFactory = Dune::Dorie::AdaptivityHandlerFactory<Traits,GFS,Parameters,FlowBoundary,U>;
 
 	Dune::MPIHelper& helper;
 	std::shared_ptr<Grid> grid;

dune/dorie/solver/operator_error_indicator.hh

@@ -13,26 +13,13 @@
 namespace Dune {
   namespace Dorie {
 
-    //! Weighting at grid interface for error estimation
-    struct EstimatorWeights
-    {
-      enum Type { weightsOn, //!< Use saturated conductivities at interface for weighting (default)
-                  weightsOff //!< No weighting
-                };
-    };
-
-    enum Impl { old, corrected, renew };
-
     /// Local operator for residual-based error estimation.
     /*
      * A call to residual() of a grid operator space will assemble
-     * the quantity \eta_T^2 for each cell. Note that the squares
-     * of the cell indicator \eta_T is stored. To compute the global
-     * error estimate sum up all values and take the square root.
+     * the quantity \eta_T for each cell.
      * Skeleton integral is mostly the same as in the DG operator.
-     *
      */
-    template<typename Traits, typename Parameter, typename FEM>
+    template<typename Traits, typename Parameter, typename Boundary>
     class FluxErrorEstimator
       : public Dune::PDELab::LocalOperatorDefaultFlags
     {
@@ -51,25 +38,27 @@ namespace Dune {
       typedef typename Traits::Domain             Domain;
       typedef typename Traits::IntersectionDomain ID;
 
-      typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
-      typedef Dune::PDELab::LocalBasisCache<LocalBasisType> Cache;
-
       // pattern assembly flags
       enum { doPatternVolume = false };
       enum { doPatternSkeleton = false };
 
       // residual assembly flags
       enum { doAlphaSkeleton  = true };
+      enum { doAlphaBoundary  = true };
 
-      Parameter& param;
+      const Parameter& param;
+      const Boundary& boundary;
       const RF penalty_factor;
       const int intorderadd;
       const int quadrature_factor;
+      RF time;
 
-      FluxErrorEstimator (const Dune::ParameterTree& config, Parameter& param_,
+      FluxErrorEstimator (const Dune::ParameterTree& config,
+        const Parameter& param_, const Boundary& boundary_,
         int intorderadd_ = 2, int quadrature_factor_ = 2)
-      : param(param_), penalty_factor(config.get<RF>("dg.penaltyFactor")),
-      intorderadd(intorderadd_), quadrature_factor(quadrature_factor_)
+      : param(param_), boundary(boundary_),
+        penalty_factor(config.get<RF>("dg.penaltyFactor")),
+        intorderadd(intorderadd_), quadrature_factor(quadrature_factor_)
       { }
 
       /// skeleton integral depending on test and ansatz functions
@@ -210,8 +199,110 @@ namespace Dune {
         r_n.accumulate(lfsv_n,0, std::sqrt( h_T_n * C_F_T_n * sum ));
       }
 
+    template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
+    void alpha_boundary (const IG& ig,
+        const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
+        R& r_s) const
+    {
+      // get polynomial degree
+      const int degree = lfsu_s.finiteElement().localBasis().order();
+      const int intorder = intorderadd + quadrature_factor * degree;
+
+      // geometric factor of penalty
+      const RF h_F = ig.inside().geometry().volume() / ig.geometry().volume(); // Houston!
+
+      // get element geometry
+      auto gtface = ig.geometry();
+
+      // container for sum of errors
+      RF sum(0.0);
+
+      // loop over quadrature points and integrate normal flux
+      for (const auto& it : quadratureRule(gtface,intorder))
+      {
+        // check boundary condition type
+        const auto bc = boundary.bc(ig.intersection(), it.position(), time);
+        if (!BoundaryCondition::isDirichlet(bc)){
+          continue;
+        }
+
+        // position of quadrature point in local coordinates of elements
+        const Domain local_s = ig.geometryInInside().global(it.position());
+
+        // retrieve unit normal vector
+        const Dune::FieldVector<DF,dim> n_F = ig.unitOuterNormal(it.position());
+
+        // evaluate basis functions
+        std::vector<Scalar> phi_s(lfsu_s.size());
+        lfsu_s.finiteElement().localBasis().evaluateFunction(local_s,phi_s);
+
+        // evaluate u
+        RF u_s = 0.;
+        for (unsigned int i = 0; i<lfsu_s.size(); i++)
+          u_s += x_s(lfsu_s,i) * phi_s[i];
+
+        // evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
+        typedef typename LFSU::Traits::FiniteElementType::Traits
+          ::LocalBasisType::Traits::JacobianType JacobianType;
+        std::vector<JacobianType> gradphi_s(lfsu_s.size());
+        lfsu_s.finiteElement().localBasis().evaluateJacobian(local_s,gradphi_s);
+
+        // transform gradients to real element
+        Dune::FieldMatrix<DF,dim,dim> jac;
+        std::vector<Vector> tgradphi_s(lfsu_s.size());
+        jac = ig.inside().geometry().jacobianInverseTransposed(local_s);
+        for (unsigned int i = 0; i<lfsu_s.size(); i++)
+          jac.mv(gradphi_s[i][0],tgradphi_s[i]);
+
+        // compute gradient of u
+        Vector gradu_s(0.);
+        for (unsigned int i = 0; i<lfsu_s.size(); i++)
+          gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
+
+        // offset for parameter queries
+        const RF eps = 1e-9;
+
+        // parameter queries: inside
+        Domain xGlobal_s = ig.inside().geometry().global(local_s);
+        xGlobal_s = xGlobal_s.axpy(-eps,n_F);
+
+        const RF head_s = param.headMillerToRef(u_s,xGlobal_s);
+        const RF satCond_s = param.condRefToMiller(param.K(xGlobal_s),xGlobal_s);
+        const RF saturation_s = param.saturation(head_s,xGlobal_s);
+        const RF condFactor_s = param.condFactor(saturation_s,xGlobal_s);
+
+        // boundary condition at position
+        const RF g = boundary.g(ig.intersection(), it.position(), time);
+
+        // compute jump in solution
+        const RF h_jump = u_s - g;
+
+        // apply gravity vector
+        gradu_s -= param.gravity();
+
+        // penalty term
+        const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
+
+        // integration factor
+        const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
+
+        const RF head_jump = satCond_s * condFactor_s * penalty * h_jump;
+
+        sum += head_jump * head_jump * factor;
+      }
+
+      DF h_T_s = diameter(ig.inside().geometry());
+
+      DF C_P_T = 1.0 / M_PI;
+      DF C_F_T_s = h_T_s * ig.geometry().volume() / ig.inside().geometry().volume();
+      C_F_T_s *= C_P_T * (2.0 / dim + C_P_T);
+
+      // accumulate indicator
+      r_s.accumulate(lfsv_s,0, std::sqrt( h_T_s * C_F_T_s * sum ));
+    }
+
     private:
-      // probably redundant
+      //! compute diameter of a grid cell
       template<class GEO>
       typename GEO::ctype diameter (const GEO& geo) const
       {
@@ -231,6 +322,13 @@ namespace Dune {
         return hmax;
       }
 
+    public:
+      //! set operator time
+      void setTime (RF t)
+      {
+        time = t;
+      }
+
     };
 
   }