operator_error_indicator.hh 11.2 KB
Newer Older
1 2 3
#ifndef DUNE_DORIE_ERROR_INDICATOR_HH
#define DUNE_DORIE_ERROR_INDICATOR_HH

4 5 6 7 8
#include <vector>

#include <dune/common/fvector.hh>
#include <dune/common/parametertree.hh>

Santiago Ospina's avatar
Santiago Ospina committed
9
#include <dune/geometry/referenceelements.hh>
Lukas Riedel's avatar
Lukas Riedel committed
10

Santiago Ospina's avatar
Santiago Ospina committed
11 12 13 14 15 16 17
#include <dune/pdelab/common/quadraturerules.hh>
#include <dune/pdelab/localoperator/pattern.hh>
#include <dune/pdelab/localoperator/flags.hh>
#include <dune/pdelab/localoperator/idefault.hh>
#include <dune/pdelab/localoperator/defaultimp.hh>
#include <dune/pdelab/finiteelement/localbasiscache.hh>

18 19
#include <dune/dorie/common/typedefs.hh>
#include <dune/dorie/model/richards/local_operator_cfg.hh>
Dion Haefner's avatar
Dion Haefner committed
20 21 22

namespace Dune {
  namespace Dorie {
23
  namespace Operator { 
Dion Haefner's avatar
Dion Haefner committed
24 25 26 27

    /// Local operator for residual-based error estimation.
    /*
     * A call to residual() of a grid operator space will assemble
28
     * the quantity \eta_T for each cell.
Dion Haefner's avatar
Dion Haefner committed
29 30
     * Skeleton integral is mostly the same as in the DG operator.
     */
31
    template<typename Traits, typename Parameter, typename Boundary>
Dion Haefner's avatar
Dion Haefner committed
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
    class FluxErrorEstimator
      : public Dune::PDELab::LocalOperatorDefaultFlags
    {

    public:

      //typedef typename Parameter::Traits::GridTraits GridTraits;

      enum { dim = Traits::dim };

      typedef typename Traits::GridView           GridView;
      typedef typename Traits::RangeField         RF;
      typedef typename Traits::Scalar             Scalar;
      typedef typename Traits::Vector             Vector;
      typedef typename Traits::DomainField        DF;
      typedef typename Traits::Domain             Domain;
      typedef typename Traits::IntersectionDomain ID;

      // pattern assembly flags
      enum { doPatternVolume = false };
      enum { doPatternSkeleton = false };

      // residual assembly flags
      enum { doAlphaSkeleton  = true };
56
      enum { doAlphaBoundary  = true };
Dion Haefner's avatar
Dion Haefner committed
57

58 59
      const Parameter& param;
      const Boundary& boundary;
60 61
      const int intorderadd;
      const int quadrature_factor;
62
      RF penalty_factor;
63
      RF time;
Dion Haefner's avatar
Dion Haefner committed
64

65 66
      FluxErrorEstimator (const Dune::ParameterTree& config,
        const Parameter& param_, const Boundary& boundary_,
67
        RichardsDGMethod method_ = RichardsDGMethod::SIPG,
Dion Haefner's avatar
Dion Haefner committed
68
        int intorderadd_ = 2, int quadrature_factor_ = 2)
69
      : param(param_), boundary(boundary_),
70
        intorderadd(intorderadd_), quadrature_factor(quadrature_factor_),
71
        penalty_factor(config.get<RF>("numerics.penaltyFactor")),
72 73
        time(0.0)
      {
74
        if (method_ == RichardsDGMethod::OBB)
75 76
          penalty_factor = 0.0;
      }
Dion Haefner's avatar
Dion Haefner committed
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99

      /// skeleton integral depending on test and ansatz functions
      /** each face is only visited ONCE!
       */
      template<typename IG,
               typename LFSU,
               typename X,
               typename LFSV,
               typename R>
      void alpha_skeleton (const IG& ig,
                           const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
                           const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
                           R& r_s, R& r_n) const
      {
        // get polynomial degree
        const int order_s  = lfsu_s.finiteElement().localBasis().order();
        const int order_n  = lfsu_n.finiteElement().localBasis().order();
        const int degree   = std::max( order_s, order_n );
        const int intorder = intorderadd + quadrature_factor * degree;

        // geometric factor of penalty
        const RF h_F = std::min(ig.inside().geometry().volume(),ig.outside().geometry().volume())/ig.geometry().volume(); // Houston!

100 101
        // get element geometry
        auto gtface = ig.geometry();
Dion Haefner's avatar
Dion Haefner committed
102 103 104 105

        // container for sum of errors
        RF sum(0.0);

106 107 108 109 110 111 112 113 114
        // bind parameterization and retrieve functions
        param.bind(ig.inside());
        const auto saturation_f_s = param.saturation_f();
        const auto conductivity_f_s = param.conductivity_f();

        param.bind(ig.outside());
        const auto saturation_f_n = param.saturation_f();
        const auto conductivity_f_n = param.conductivity_f();

Dion Haefner's avatar
Dion Haefner committed
115
        // loop over quadrature points and integrate normal flux
116
        for (const auto& it : quadratureRule(gtface,intorder))
Dion Haefner's avatar
Dion Haefner committed
117
        {
118 119
          // retrieve unit normal vector
          const Dune::FieldVector<DF,dim> n_F = ig.unitOuterNormal(it.position());
Dion Haefner's avatar
Dion Haefner committed
120 121

          // position of quadrature point in local coordinates of elements
122 123
          const Domain local_s = ig.geometryInInside().global(it.position());
          const Domain local_n = ig.geometryInOutside().global(it.position());
Dion Haefner's avatar
Dion Haefner committed
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163

          // evaluate basis functions
          std::vector<Scalar> phi_s(lfsu_s.size());
          std::vector<Scalar> phi_n(lfsu_n.size());
          lfsu_s.finiteElement().localBasis().evaluateFunction(local_s,phi_s);
          lfsu_n.finiteElement().localBasis().evaluateFunction(local_n,phi_n);

          // evaluate u
          RF u_s = 0., u_n = 0.;
          for (unsigned int i = 0; i<lfsu_s.size(); i++)
            u_s += x_s(lfsu_s,i) * phi_s[i];
          for (unsigned int i = 0; i<lfsu_n.size(); i++)
            u_n += x_n(lfsu_n,i) * phi_n[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());
          std::vector<JacobianType> gradphi_n(lfsu_n.size());
          lfsu_s.finiteElement().localBasis().evaluateJacobian(local_s,gradphi_s);
          lfsu_n.finiteElement().localBasis().evaluateJacobian(local_n,gradphi_n);

          // transform gradients to real element
          Dune::FieldMatrix<DF,dim,dim> jac;
          std::vector<Vector> tgradphi_s(lfsu_s.size());
          std::vector<Vector> tgradphi_n(lfsu_n.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]);
          jac = ig.outside().geometry().jacobianInverseTransposed(local_n);
          for (unsigned int i = 0; i<lfsu_n.size(); i++)
            jac.mv(gradphi_n[i][0],tgradphi_n[i]);

          // compute gradient of u
          Vector gradu_s(0.), gradu_n(0.);
          for (unsigned int i = 0; i<lfsu_s.size(); i++)
            gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
          for (unsigned int i = 0; i<lfsu_n.size(); i++)
            gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);

164 165 166
          // retrieve conductivity
          const RF cond_s = conductivity_f_s(saturation_f_s(u_s));
          const RF cond_n = conductivity_f_n(saturation_f_n(u_n));
Dion Haefner's avatar
Dion Haefner committed
167 168 169 170 171 172 173 174 175 176 177 178

          // compute jump in solution
          const RF h_jump = u_s - u_n;

          // apply gravity vector
          gradu_s -= param.gravity();
          gradu_n -= param.gravity();

          // penalty term
          const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);

          // integration factor
179
          const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
Dion Haefner's avatar
Dion Haefner committed
180

181 182
          const RF grad_jump = (cond_s * (n_F * gradu_s))
                               - (cond_n * (n_F * gradu_n));
183

184
          const RF gamma = 2.0 * cond_s * cond_n / ( cond_s + cond_n );
185
          const RF head_jump = gamma * penalty * h_jump;
186

187
          const RF total_jump = grad_jump/2.0 + head_jump;
188

189
          sum += total_jump * total_jump * factor;
190
        }
Dion Haefner's avatar
Dion Haefner committed
191

192 193
        DF h_T_s = diameter(ig.inside().geometry());
        DF h_T_n = diameter(ig.outside().geometry());
Dion Haefner's avatar
Dion Haefner committed
194

195
        DF C_P_T = 1.0 / M_PI;
196 197 198 199
        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);
        DF C_F_T_n = h_T_n * ig.geometry().volume() / ig.outside().geometry().volume();
        C_F_T_n *= C_P_T * (2.0 / dim + C_P_T);
Dion Haefner's avatar
Dion Haefner committed
200 201

        // accumulate indicator
202 203
        r_s.accumulate(lfsv_s,0, h_T_s * C_F_T_s * sum );
        r_n.accumulate(lfsv_n,0, h_T_n * C_F_T_n * sum );
Dion Haefner's avatar
Dion Haefner committed
204 205
      }

206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
    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);

224 225 226 227 228
      // bind parameterization and retrieve functions
      param.bind(ig.inside());
      const auto saturation_f_s = param.saturation_f();
      const auto conductivity_f_s = param.conductivity_f();

229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261
      // 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());

        // 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];

        // 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;

        // 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());

262 263 264 265
        // query conductivity
        const RF cond_s = conductivity_f_s(saturation_f_s(u_s));

        const RF head_jump = cond_s * penalty * h_jump;
266 267 268 269 270 271 272 273 274 275 276

        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
277
      r_s.accumulate(lfsv_s,0, h_T_s * C_F_T_s * sum );
278 279
    }

Dion Haefner's avatar
Dion Haefner committed
280
    private:
281
      //! compute diameter of a grid cell
Dion Haefner's avatar
Dion Haefner committed
282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
      template<class GEO>
      typename GEO::ctype diameter (const GEO& geo) const
      {
        typedef typename GEO::ctype DF;
        DF hmax = -1.0E00;
        const int dim = GEO::coorddimension;
        for (int i=0; i<geo.corners(); i++)
          {
            Dune::FieldVector<DF,dim> xi = geo.corner(i);
            for (int j=i+1; j<geo.corners(); j++)
              {
                Dune::FieldVector<DF,dim> xj = geo.corner(j);
                xj -= xi;
                hmax = std::max(hmax,xj.two_norm());
              }
          }
        return hmax;
      }

301 302 303 304 305 306 307
    public:
      //! set operator time
      void setTime (RF t)
      {
        time = t;
      }

Dion Haefner's avatar
Dion Haefner committed
308 309
    };

310 311 312
  } // namespace Operator
  } // namespace Dorie
} // namespace Dune
313

314
#endif // DUNE_DORIE_ERROR_INDICATOR_HH