local_operator_DG.hh 36.9 KB
Newer Older
Dion Haefner's avatar
Dion Haefner committed
1 2 3 4
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_DORIE_RICHARDS_OPERATOR_HH
#define DUNE_DORIE_RICHARDS_OPERATOR_HH

5
#include <vector>
6
#include <memory>
7 8 9 10

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

Santiago Ospina's avatar
Santiago Ospina committed
11
#include <dune/geometry/referenceelements.hh>
Lukas Riedel's avatar
Lukas Riedel committed
12 13 14

#include <dune/grid/common/mcmgmapper.hh>

Santiago Ospina's avatar
Santiago Ospina committed
15 16 17 18 19 20 21
#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>

22
#include <dune/dorie/common/typedefs.hh>
23
#include <dune/dorie/model/richards/local_operator_cfg.hh>
Dion Haefner's avatar
Dion Haefner committed
24 25 26 27 28 29 30 31

namespace Dune {
namespace Dorie {
namespace Operator {

//! Operator internal BC handling
struct BCType
{
32
  enum Type { Neumann, //!< Fixed flux at boundary
33
              Dirichlet, //!< Fixed matric head at boundary
Dion Haefner's avatar
Dion Haefner committed
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
              none //!< No BC specified
            };
};

/// A local operator for solving the Richards equation using upwinding.
/** Written by Ole Klein, IWR, Heidelberg and modified by the
 *  DORiE developers, IUP, Heidelberg \n
 *  In theory it is possible that one and the same local operator is
 *  called first with a finite element of one type and later with a
 *  finite element of another type.  Since finite elements of different
 *  type will usually produce different results for the same local
 *  coordinate they cannot share a cache.  Here we use a vector of caches
 *  to allow for different orders of the shape functions, which should be
 *  enough to support p-adaptivity.  (Another likely candidate would be
 *  differing geometry types, i.e. hybrid meshes.)
49
 *
50 51
 * @ingroup    LocalOperators
 * @ingroup    RichardsModel
Dion Haefner's avatar
Dion Haefner committed
52 53 54 55 56 57 58 59 60 61 62 63 64 65
 */
template<typename Traits, typename Parameter, typename Boundary, typename SourceTerm, typename FEM, bool adjoint>
  class RichardsDGSpatialOperator
  : public Dune::PDELab::NumericalJacobianVolume<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint> >,
  public Dune::PDELab::NumericalJacobianSkeleton<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint> >,
  public Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint> >,
  public Dune::PDELab::FullSkeletonPattern,
  public Dune::PDELab::FullVolumePattern,
  public Dune::PDELab::LocalOperatorDefaultFlags,
  public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<typename Traits::RangeField>
{
protected:

  enum { dim = Traits::dim };
66
  static_assert(dim!=1, "RichardsDGSpatialOperator does not work for world dimension 1!");
67
  
Dion Haefner's avatar
Dion Haefner committed
68 69 70 71 72 73 74 75 76 77 78 79 80
  typedef typename Traits::RangeField         RF;
  typedef typename Traits::Scalar             Scalar;
  typedef typename Traits::Vector             Vector;
  typedef typename Traits::Tensor             Tensor;
  typedef typename Traits::DomainField        DF;
  typedef typename Traits::Domain             Domain;
  typedef typename Traits::IntersectionDomain ID;
  typedef typename Traits::Element            Element;


  typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
  typedef Dune::PDELab::LocalBasisCache<LocalBasisType> Cache;

81 82 83
  std::shared_ptr<const Parameter> param; //!< class providing equation parameters
  std::shared_ptr<const Boundary> boundary; //!< class providing boundary conditions
  std::shared_ptr<const SourceTerm> sourceTerm; //!< class providing source term information
84 85 86
  const RichardsDGMethod method; //!< Interior penalty type
  const RichardsUpwinding upwinding; //!< DG method weights switch
  const RichardsDGWeights weights; //!< Gradient weighting
Dion Haefner's avatar
Dion Haefner committed
87 88 89

  const int intorderadd; //!< integration order addend
  const int quadrature_factor; //!< factor to FEM integration order
90 91 92
  RF penalty_factor; //!< Interior penalty factor
  RF theta; //!< Symmetry factor
  RF time;  //!< operator internal time for BC, parameter, and source term queries
Dion Haefner's avatar
Dion Haefner committed
93 94 95

  std::vector<Cache> cache; //!< Local basis cache

96 97
  enum class LOPCase {DG,RTVolume,RTSkeleton}; //!< Local operator case

Dion Haefner's avatar
Dion Haefner committed
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114
public:

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

  // residual assembly flags
  enum { doAlphaVolume  = true };
  enum { doAlphaSkeleton  = true };
  enum { doAlphaBoundary  = true };
  enum { doLambdaVolume   = true };

  /// Create operator.
  /**
   *  \param param_ Object providing equation parameters
   *  \param boundary_ Object prividing boundary conditions
   *  \param sourceTerm_ Object providing source term information
115 116 117
   *  \param method_ DG method type \see RichardsDGMethod
   *  \param upwinding_ Upwinding type \see RichardsUpwinding
   *  \param weights_ DG weigths switch \see RichardsDGWeights
Dion Haefner's avatar
Dion Haefner committed
118 119 120
   *  \param intorderadd_ Addend to integration order
   *  \param quadrature_factor_ Factor to FEM integration order
   */
121 122
  RichardsDGSpatialOperator (
    const Dune::ParameterTree& config,
123 124 125
    std::shared_ptr<const Parameter> param_,
    std::shared_ptr<const Boundary> boundary_,
    std::shared_ptr<const SourceTerm> sourceTerm_,
126 127 128
    const RichardsDGMethod method_ = RichardsDGMethod::SIPG,
    const RichardsUpwinding upwinding_ = RichardsUpwinding::none,
    const RichardsDGWeights weights_ = RichardsDGWeights::weightsOn,
129 130 131 132
    int intorderadd_ = 2,
    int quadrature_factor_ = 2
  ) :
    Dune::PDELab::NumericalJacobianVolume<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
Dion Haefner's avatar
Dion Haefner committed
133 134
    Dune::PDELab::NumericalJacobianSkeleton<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
    Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
135 136 137 138 139 140 141 142
    param(param_),
    boundary(boundary_),
    sourceTerm(sourceTerm_),
    method(method_),
    upwinding(upwinding_),
    weights(weights_),
    intorderadd(intorderadd_),
    quadrature_factor(quadrature_factor_),
143
    penalty_factor(config.get<RF>("numerics.penaltyFactor")),
144
    time(0.0),
Dion Haefner's avatar
Dion Haefner committed
145
    cache(20)
Dion Haefner's avatar
Dion Haefner committed
146
  {
147
    theta = 0.; // IIP
148
    if (method == RichardsDGMethod::NIPG || method == RichardsDGMethod::OBB)
149 150 151
      theta = 1.;
    else if (method == RichardsDGMethod::SIPG)
      theta = -1.;
152

153
    if (method == RichardsDGMethod::OBB)
154
      penalty_factor = 0.;
Dion Haefner's avatar
Dion Haefner committed
155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
  }

  /// Volume integral depending on test and ansatz functions
  /** Richards Equation - Strong formulation of volume integral
   *  \f[
   *    \int_{\Omega_e} K (\Theta (h_m^*)) \cdot (\nabla u - \mathbf{\hat{g}}) \cdot \nabla v \; dx
   *  \f]
   *  Miller Similarity Transformations:
   *  \f{eqnarray*}{
   *    h_m \; &\xrightarrow{\text{MS}^{-1}} \; & h_m^* \\
   *    K^*(\Theta(h_m^*)) \; &\xrightarrow{\text{MS}} \;  & K(\Theta(h_m^*))
   *  \f}
   */
  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
    void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
    {
171
      // Switches between local and global interface in test local space
172 173
      using FETest = typename LFSV::Traits::FiniteElementType;
      using BasisTest = typename FETest::Traits::LocalBasisType;
174 175 176 177 178 179 180

      /* If range dimension of the test local function space is 1 we assume that
       * it is a P_k/Q_k space (probably conforming), if it is equal to world 
       * dimension we assume that it is a vector function space representing a 
       * gradient (e.g. Raviart Thomas) which is non conforming with respect to 
       * the ansatz space.*/

181
      constexpr int dimRangeTest = BasisTest::Traits::dimRange;
182
      constexpr LOPCase lopcase = (dimRangeTest==dim) ? LOPCase::RTVolume : LOPCase::DG;
183

Dion Haefner's avatar
Dion Haefner committed
184 185 186 187 188 189
      const int order = lfsu.finiteElement().localBasis().order();
      const int intorder = intorderadd + quadrature_factor * order;

      // get element geomentry
      auto gt = eg.geometry();

190
      // bind parameterization and retrieve functions
191 192 193
      param->bind(eg.entity());
      const auto saturation_f = param->saturation_f();
      const auto conductivity_f = param->conductivity_f();
194

Dion Haefner's avatar
Dion Haefner committed
195 196 197
      // loop over quadrature points
      for (const auto& it : quadratureRule(gt,intorder))
      {
198
        // evaluate position in element local and global coordinates
199
        const auto p_local = it.position();
200

Dion Haefner's avatar
Dion Haefner committed
201
        // evaluate basis functions
202
        auto& phiu = cache[order].evaluateFunction(p_local,lfsu.finiteElement().localBasis());
Dion Haefner's avatar
Dion Haefner committed
203 204 205 206

        // evaluate u
        RF u = 0.;
        for (unsigned int i = 0; i<lfsu.size(); i++)
207
          u += x(lfsu,i) * phiu[i];
Dion Haefner's avatar
Dion Haefner committed
208

209
        // evaluate gradient of basis function
210
        const auto& js = cache[order].evaluateJacobian(p_local,lfsu.finiteElement().localBasis());
Dion Haefner's avatar
Dion Haefner committed
211 212

        // transform gradients to real element
213
        const auto jac = gt.jacobianInverseTransposed(p_local);
214
        std::vector<Vector> gradphiu(lfsu.size());
Dion Haefner's avatar
Dion Haefner committed
215
        for (unsigned int i = 0; i<lfsu.size(); i++)
216 217
          jac.mv(js[i][0],gradphiu[i]);

218
        std::vector<Vector> gradphiv(lfsv.size());
219

220
        if constexpr (lopcase == LOPCase::RTVolume) {
221 222
          // (we assume lfsv = RTk)
          lfsv.finiteElement().localBasis().evaluateFunction(p_local,gradphiv);
223
          // Piola transformation
224 225 226 227
          auto BTransposed = gt.jacobianTransposed(p_local);
          auto detB = BTransposed.determinant();
          for (unsigned int i=0; i<lfsv.size(); i++) {
            Vector y(gradphiv[i]);
228
            y /= detB;
229 230
            BTransposed.mtv(y,gradphiv[i]);
          }
231
        } else if constexpr (lopcase == LOPCase::DG) {
232 233
          // (we assume Galerkin method lfsu = lfsv)
          gradphiv = gradphiu;
234 235 236
        } else {
          static_assert(Dune::AlwaysTrue<EG>::value,
            "alpha_volume() does not support RTSkeleton case!");
237
        }
Dion Haefner's avatar
Dion Haefner committed
238 239 240 241

        // compute gradient of u
        Vector gradu(0.);
        for (unsigned int i = 0; i<lfsu.size(); i++)
242
          gradu.axpy(x(lfsu,i),gradphiu[i]);
Dion Haefner's avatar
Dion Haefner committed
243

244
        gradu -= param->gravity();
Dion Haefner's avatar
Dion Haefner committed
245

246 247
        // compute conductivity
        const RF cond = conductivity_f(saturation_f(u));
Dion Haefner's avatar
Dion Haefner committed
248 249

        // integration factor
250
        const RF factor = it.weight() * gt.integrationElement(p_local);
Dion Haefner's avatar
Dion Haefner committed
251 252

        // update residual
253 254 255 256 257
        if constexpr (lopcase != LOPCase::RTSkeleton) {
          
          // discrete gradient sign: Negative for lifting
          double dg_sign = (lopcase == LOPCase::RTVolume) ? -1. : 1.;
          
258 259
          for (unsigned int i = 0; i<lfsv.size(); i++)
            r.accumulate(lfsv,i, dg_sign * cond * (gradu*gradphiv[i]) * factor);
260
        }
Dion Haefner's avatar
Dion Haefner committed
261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289
      }
    }

  /// Skeleton integral depending on test and ansatz functions
  /** Each face is only visited ONCE!\n
   * Weights for symmetry term:
   *  \f[
   *    \omega_x = \frac{ K_{0,y} }{ K_{0,x} + K_{0,y}} \; , \quad x,y \in \lbrace s,n \rbrace , \; x \neq y
   *  \f]
   * Solution jump over face:
   *  \f[
   *    \Delta u = u_s - u_n
   *  \f]
   * Penalty factor:
   *  \f{eqnarray*}{
   *    p &= \mathsf{p} \cdot k (k+d-1) / h_f \\
   *    h_f &= \frac{\min \lbrace \Omega_s,\Omega_n \rbrace }{\Omega_i}\\
   *    d \; &: \text{spatial dimension}\\
   *    k \; &: \text{polynomial order}\\
   *    \Omega_s,\Omega_n \in \mathbb{R}^d \; &: \text{element volumes}\\
   *    \Omega_i \in \mathbb{R}^{d-1} \; &: \text{interface surface}
   *  \f}
   */
  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
    {
290
      // Switches between local and global interface in test local space
291 292
      using FETest = typename LFSV::Traits::FiniteElementType;
      using BasisTest = typename FETest::Traits::LocalBasisType;
293 294 295 296 297 298

      /* If domain dimension of the test local function space is dim we assume
       * that it is a P_k/Q_k space (probably conforming), if it is dim-1 we
       * assume that it is a skeleton finite element function space
       * which is non conforming with respect to the ansatz space.*/

299 300
      constexpr int dimDomainLocalTest = BasisTest::Traits::dimDomain;
      constexpr int dimRangeTest = BasisTest::Traits::dimRange;
301 302
      constexpr LOPCase lopcase = (dimDomainLocalTest==dim-1) ? LOPCase::RTSkeleton : 
                                  (dimRangeTest==dim ? LOPCase::RTVolume : LOPCase::DG);
303

Dion Haefner's avatar
Dion Haefner committed
304 305 306 307 308 309 310 311 312
      // 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;

      // get element geometry
      auto gtface = ig.geometry();

313
      // bind parameterization and retrieve functions
314 315 316
      param->bind(ig.inside());
      const auto saturation_f_s = param->saturation_f();
      const auto conductivity_f_s = param->conductivity_f();
317 318

      // cache saturated conductivity
319
      const RF satCond_s = conductivity_f_s(1.0);
320

321 322 323
      param->bind(ig.outside());
      const auto saturation_f_n = param->saturation_f();
      const auto conductivity_f_n = param->conductivity_f();
324

325
      const RF satCond_n = conductivity_f_n(1.0);
326

Dion Haefner's avatar
Dion Haefner committed
327 328 329 330 331 332 333
      // geometric factor of penalty
      const RF h_F = std::min(ig.inside().geometry().volume(),ig.outside().geometry().volume())/ig.geometry().volume(); // Houston!

      // loop over quadrature points and integrate normal flux
      for (const auto& it : quadratureRule(gtface,intorder))
      {
        // position of quadrature point in local coordinates of elements
334 335
        const auto p_local_s = ig.geometryInInside().global(it.position());
        const auto p_local_n = ig.geometryInOutside().global(it.position());
Dion Haefner's avatar
Dion Haefner committed
336

337 338 339
        // face normal vector
        const Domain normal = ig.unitOuterNormal(it.position());

340 341 342
        // evaluate basis functions of trial function
        const auto& phiu_s = cache[order_s].evaluateFunction(p_local_s,lfsu_s.finiteElement().localBasis());
        const auto& phiu_n = cache[order_n].evaluateFunction(p_local_n,lfsu_n.finiteElement().localBasis());
Dion Haefner's avatar
Dion Haefner committed
343

344 345
        std::vector<Scalar> phiv_s(lfsv_s.size());
        std::vector<Scalar> phiv_n(lfsv_n.size());
Dion Haefner's avatar
Dion Haefner committed
346

347 348 349
        // evaluate gradient of basis function
        const auto& jsu_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
        const auto& jsu_n = cache[order_n].evaluateJacobian(p_local_n,lfsu_n.finiteElement().localBasis());
Dion Haefner's avatar
Dion Haefner committed
350 351 352

        // transform gradients to real element
        Tensor jac;
353 354
        std::vector<Vector> gradphiu_s(lfsu_s.size());
        std::vector<Vector> gradphiu_n(lfsu_n.size());
355
        jac = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
Dion Haefner's avatar
Dion Haefner committed
356
        for (unsigned int i = 0; i<lfsu_s.size(); i++)
357
          jac.mv(jsu_s[i][0],gradphiu_s[i]);
358
        jac = ig.outside().geometry().jacobianInverseTransposed(p_local_n);
Dion Haefner's avatar
Dion Haefner committed
359
        for (unsigned int i = 0; i<lfsu_n.size(); i++)
360 361 362 363 364
          jac.mv(jsu_n[i][0],gradphiu_n[i]);

        std::vector<Vector> gradphiv_s(lfsv_s.size());
        std::vector<Vector> gradphiv_n(lfsv_n.size());

365
        // evaluate basis functions of test function
366
        if constexpr (lopcase == LOPCase::RTSkeleton) {
367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383
          // Evaluate position in local coordinates of the inside entity of codim 1
          using LocalCoordinate = typename IG::LocalGeometry::LocalCoordinate;
          LocalCoordinate sub_p_local_s, sub_p_local_n;

          if (ig.intersection().conforming()) {
            sub_p_local_s = it.position();
            sub_p_local_n = it.position();
          } else {
            auto sub_entity_s_id = ig.intersection().indexInInside();
            auto sub_entity_s = ig.inside().template subEntity<1>(sub_entity_s_id);
            sub_p_local_s = sub_entity_s.geometry().local(gtface.global(it.position()));

            auto sub_entity_n_id = ig.intersection().indexInOutside();
            auto sub_entity_n = ig.outside().template subEntity<1>(sub_entity_n_id);
            sub_p_local_n = sub_entity_n.geometry().local(gtface.global(it.position()));
          }

384
          // (we assume non-Galerkin method lfsu != lfsv)
385 386
          lfsv_s.finiteElement().localBasis().evaluateFunction(sub_p_local_s,phiv_s);
          lfsv_n.finiteElement().localBasis().evaluateFunction(sub_p_local_n,phiv_n);
387
        } else if constexpr (lopcase == LOPCase::RTVolume) {
388 389 390 391
          // (we assume lfsv = RTk)
          lfsv_s.finiteElement().localBasis().evaluateFunction(p_local_s,gradphiv_s);
          lfsv_n.finiteElement().localBasis().evaluateFunction(p_local_n,gradphiv_n);
          // Piola transformation
392 393
          auto BTransposed_s = ig.inside().geometry().jacobianTransposed(p_local_s);
          auto detB_s = BTransposed_s.determinant();
394 395
          for (unsigned int i=0; i<lfsv_s.size(); i++) {
            Vector y_s(gradphiv_s[i]);
396 397
            y_s /= detB_s;
            BTransposed_s.mtv(y_s,gradphiv_s[i]);
398
          }
399
          auto BTransposed_n = ig.outside().geometry().jacobianTransposed(p_local_n);
400
          auto detB_n = BTransposed_n.determinant();
401 402
          for (unsigned int i=0; i<lfsv_n.size(); i++) {
            Vector y_n(gradphiv_n[i]);
403 404
            y_n /= detB_n;
            BTransposed_n.mtv(y_n,gradphiv_n[i]);
405 406
          }

407 408 409 410
        } else {
          // (we assume Galerkin method lfsu = lfsv)
          phiv_s = phiu_s;
          phiv_n = phiu_n;
411 412
          gradphiv_s = gradphiu_s;
          gradphiv_n = gradphiu_n;
413
        }
Dion Haefner's avatar
Dion Haefner committed
414 415 416 417

        // evaluate u
        RF u_s = 0., u_n = 0.;
        for (unsigned int i = 0; i<lfsu_s.size(); i++)
418
          u_s += x_s(lfsu_s,i) * phiu_s[i];
Dion Haefner's avatar
Dion Haefner committed
419
        for (unsigned int i = 0; i<lfsu_n.size(); i++)
420
          u_n += x_n(lfsu_n,i) * phiu_n[i];
Dion Haefner's avatar
Dion Haefner committed
421 422 423 424

        // compute gradient of u
        Vector gradu_s(0.), gradu_n(0.);
        for (unsigned int i = 0; i<lfsu_s.size(); i++)
425
          gradu_s.axpy(x_s(lfsu_s,i),gradphiu_s[i]);
Dion Haefner's avatar
Dion Haefner committed
426
        for (unsigned int i = 0; i<lfsu_n.size(); i++)
427
          gradu_n.axpy(x_n(lfsu_n,i),gradphiu_n[i]);
Dion Haefner's avatar
Dion Haefner committed
428 429

        // update with gravity vector
430 431
        gradu_s -= param->gravity();
        gradu_n -= param->gravity();
Dion Haefner's avatar
Dion Haefner committed
432 433 434 435

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

436 437 438
        // conductivities
        const RF cond_s = conductivity_f_s(saturation_f_s(u_s));
        const RF cond_n = conductivity_f_n(saturation_f_n(u_n));
439

Dion Haefner's avatar
Dion Haefner committed
440 441
        // compute weights
        RF omega_s = 0.5, omega_n = 0.5;
442
        if (weights == RichardsDGWeights::weightsOn)
443
        {
444
          if (upwinding == RichardsUpwinding::none)
445
          {
446 447
            omega_s = cond_n / ((cond_s) + (cond_n));
            omega_n = cond_s / ((cond_s) + (cond_n));
448
          }
449 450
          else if (upwinding == RichardsUpwinding::semiUpwind
            || upwinding == RichardsUpwinding::fullUpwind)
451 452 453 454
          {
            omega_s = satCond_n / (satCond_s + satCond_n);
            omega_n = satCond_s / (satCond_s + satCond_n);
          }
Dion Haefner's avatar
Dion Haefner committed
455 456 457 458 459
        }

        // integration factor
        const RF factor = it.weight() * ig.geometry().integrationElement(it.position());

460
        // penalty factor
Dion Haefner's avatar
Dion Haefner committed
461 462
        const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);

463 464 465 466
        // relative conductivites
        const RF relCond_s = cond_s / satCond_s;
        const RF relCond_n = cond_n / satCond_n;

Dion Haefner's avatar
Dion Haefner committed
467
        // compute numerical flux
468 469
        RF numFlux_s(0.);
        RF numFlux_n(0.);
470 471 472 473
        
        numFlux_s = numFlux_n = penalty * jump;
        numFlux_s -= gradu_s * normal;
        numFlux_n -= gradu_n * normal;
474 475 476
        numFlux_s *= satCond_s;
        numFlux_n *= satCond_n;

477
        if (upwinding == RichardsUpwinding::none)
478 479 480 481 482
        {
          numFlux_s *= relCond_s;
          numFlux_n *= relCond_n;
        }
        const RF numFlux = omega_s * numFlux_s + omega_n * numFlux_n;
Dion Haefner's avatar
Dion Haefner committed
483

484 485
        // Upwind in relative conductivity
        RF relCond_upwind=0.0, satCond_upwind=0.0;
486 487
        if (upwinding == RichardsUpwinding::semiUpwind
          || upwinding == RichardsUpwinding::fullUpwind)
488
        {
489
          RF cond_upwind;
490 491
          if (numFlux > 0)
          {
492
            cond_upwind = conductivity_f_s(saturation_f_s(u_s));
493
            satCond_upwind = satCond_s;
494 495 496
          }
          else
          {
497
            cond_upwind = conductivity_f_n(saturation_f_n(u_n));
498 499
            satCond_upwind = satCond_n;
          }
500 501

          relCond_upwind = cond_upwind / satCond_upwind;
502
        }
Dion Haefner's avatar
Dion Haefner committed
503

504
        if constexpr (lopcase != LOPCase::RTVolume)
505
        {
506 507 508 509
          // update residual (flux term)
          // diffusion term
          // consistency term
          // + penalty term
510
          if (upwinding == RichardsUpwinding::none)
511 512
          {
            for (unsigned int i = 0; i < lfsv_s.size(); i++)
513
              r_s.accumulate(lfsv_s, i,   numFlux * phiv_s[i] * factor);
514
            for (unsigned int i = 0; i < lfsv_n.size(); i++)
515
              r_n.accumulate(lfsv_n, i, - numFlux * phiv_n[i] * factor);
516
          }
517
          else if (upwinding == RichardsUpwinding::semiUpwind)
518 519
          {
            for (unsigned int i = 0; i<lfsv_s.size(); i++)
520
              r_s.accumulate(lfsv_s,i,   relCond_upwind * numFlux * phiv_s[i] * factor);
521
            for (unsigned int i = 0; i<lfsv_n.size(); i++)
522
              r_n.accumulate(lfsv_n,i, - relCond_upwind * numFlux * phiv_n[i] * factor);
523
          }
524
          else if (upwinding == RichardsUpwinding::fullUpwind)
525
          {
526 527 528
            // Upwind numerical flux
            double numFlux_upwind = (numFlux > 0) ? numFlux_s : numFlux_n;

529
            for (unsigned int i = 0; i < lfsv_s.size(); i++)
530
              r_s.accumulate(lfsv_s, i,   relCond_upwind * numFlux_upwind * phiv_s[i] * factor);
531
            for (unsigned int i = 0; i < lfsv_n.size(); i++)
532
              r_n.accumulate(lfsv_n, i, - relCond_upwind * numFlux_upwind * phiv_n[i] * factor);
533
          }
534
        }
Dion Haefner's avatar
Dion Haefner committed
535

536
        if constexpr (lopcase != LOPCase::RTSkeleton)
537
        {
538 539 540
          // update residual (symmetry term)
          // (non-)symmetric IP term
          // symmetry term
541 542 543 544
          
          // discrete gradient sign: Negative for lifting
          double dg_sign = (lopcase == LOPCase::RTVolume) ? -1. : 1.;

545
          if (upwinding == RichardsUpwinding::none)
546 547 548 549 550 551
          {
            for (unsigned int i = 0; i < lfsv_s.size(); i++)
              r_s.accumulate(lfsv_s, i, dg_sign * theta * omega_s * (gradphiv_s[i] * normal) * relCond_s * satCond_s * jump * factor);
            for (unsigned int i = 0; i < lfsv_n.size(); i++)
              r_n.accumulate(lfsv_n, i, dg_sign * theta * omega_n * (gradphiv_n[i] * normal) * relCond_n * satCond_n * jump * factor);
          }
552
          else if (upwinding == RichardsUpwinding::semiUpwind)
553 554 555 556 557 558
          {
            for (unsigned int i = 0; i<lfsv_s.size(); i++)
              r_s.accumulate(lfsv_s,i, dg_sign * theta * omega_s * (gradphiv_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
            for (unsigned int i = 0; i<lfsv_n.size(); i++)
              r_n.accumulate(lfsv_n,i, dg_sign * theta * omega_n * (gradphiv_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
          }
559
          else if (upwinding == RichardsUpwinding::fullUpwind)
560 561 562 563 564 565
          {
            for (unsigned int i = 0; i < lfsv_s.size(); i++)
              r_s.accumulate(lfsv_s, i, dg_sign * theta * (gradphiv_s[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
            for (unsigned int i = 0; i < lfsv_n.size(); i++)
              r_n.accumulate(lfsv_n, i, dg_sign * theta * (gradphiv_n[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
          }
566
        }
Dion Haefner's avatar
Dion Haefner committed
567 568 569 570 571 572 573 574 575 576 577 578
      }
    }

  /**
   * @brief Boundary integral depending on test and ansatz functions
   * \todo there is additional code for the limited influx BC, which is currently commented out
   */
  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
    {
579
      // Switches between local and global interface in test local space
580 581
      using FETest = typename LFSV::Traits::FiniteElementType;
      using BasisTest = typename FETest::Traits::LocalBasisType;
582 583 584 585 586 587

      /* If domain dimension of the test local function space is dim we assume
       * that it is a P_k/Q_k space (probably conforming), if it is dim-1 we
       * assume that it is a skeleton finite element function space
       * which is non conforming with respect to the ansatz space.*/

588 589
      constexpr int dimDomainLocalTest = BasisTest::Traits::dimDomain;
      constexpr int dimRangeTest = BasisTest::Traits::dimRange;
590 591
      constexpr LOPCase lopcase = (dimDomainLocalTest==dim-1) ? LOPCase::RTSkeleton : 
                                  (dimRangeTest==dim ? LOPCase::RTVolume : LOPCase::DG);
592

Dion Haefner's avatar
Dion Haefner committed
593 594 595 596 597 598 599 600 601 602 603
      // get polynomial degree
      const int order_s  = lfsu_s.finiteElement().localBasis().order();
      const int degree   = order_s;
      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();

604
      // bind parameterization and retrieve functions
605 606 607
      param->bind(ig.inside());
      const auto saturation_f_s = param->saturation_f();
      const auto conductivity_f_s = param->conductivity_f();
608 609

      // cache saturated conductivity
610
      const RF satCond_s = conductivity_f_s(1.0);
611

Dion Haefner's avatar
Dion Haefner committed
612 613 614 615
      // loop over quadrature points and integrate normal flux
      for (const auto& it : quadratureRule(gtface,intorder))
      {
        // position of quadrature point in local coordinates of elements
616
        const auto p_local_s = ig.geometryInInside().global(it.position());
Dion Haefner's avatar
Dion Haefner committed
617

618 619
        // evaluate basis functions of trial function
        const auto& phiu_s = cache[order_s].evaluateFunction(p_local_s,lfsu_s.finiteElement().localBasis());
Dion Haefner's avatar
Dion Haefner committed
620

621 622
        std::vector<Scalar> phiv_s(lfsv_s.size());

623 624
        // evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
        const auto& jsu_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
625

626 627 628 629 630
        // transform gradients to real element
        const auto js = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
        std::vector<Vector> gradphiu_s(lfsu_s.size());
        for (unsigned int i = 0; i<lfsu_s.size(); i++)
          js.mv(jsu_s[i][0],gradphiu_s[i]);
631 632 633

        std::vector<Vector> gradphiv_s(lfsv_s.size());

634
        // evaluate basis functions of test function
635
        if constexpr (lopcase == LOPCase::RTSkeleton) {
636 637
          // (we assume non-Galerkin method lfsu != lfsv)
          lfsv_s.finiteElement().localBasis().evaluateFunction(it.position(),phiv_s);
638
        } else if constexpr (lopcase == LOPCase::RTVolume) {
639 640 641
          // (we assume lfsv = RTk)
          lfsv_s.finiteElement().localBasis().evaluateFunction(p_local_s,gradphiv_s);
          // Piola transformation
642 643
          auto BTransposed_s = ig.inside().geometry().jacobianTransposed(p_local_s);
          auto detB_s = BTransposed_s.determinant();
644 645
          for (unsigned int i=0; i<lfsv_s.size(); i++) {
            Vector y_s(gradphiv_s[i]);
646 647
            y_s /= detB_s;
            BTransposed_s.mtv(y_s,gradphiv_s[i]);
648
          }
649 650 651
        } else {
          // (we assume Galerkin method lfsu = lfsv)
          phiv_s = phiu_s;
652
          gradphiv_s = gradphiu_s;
653
        }
Dion Haefner's avatar
Dion Haefner committed
654 655 656 657 658 659
        // integration factor
        const RF factor = it.weight() * ig.geometry().integrationElement(it.position());

        // evaluate u
        RF u_s = 0.;
        for (unsigned int i = 0; i<lfsu_s.size(); i++)
660
          u_s += x_s(lfsu_s,i) * phiu_s[i];
Dion Haefner's avatar
Dion Haefner committed
661

662
        // query the type of boundary condition
Dion Haefner's avatar
Dion Haefner committed
663
        // set the internal bcType flag
664
        const auto bc = boundary->bc(ig.intersection(),it.position(),time);
Dion Haefner's avatar
Dion Haefner committed
665 666 667 668 669 670
        auto bcType = BCType::none;
        RF normal_flux;

        if (BoundaryCondition::isNeumann(bc))
        {
          // evaluate flux boundary condition
671
          normal_flux = boundary->j(ig.intersection(),it.position(),time);
Dion Haefner's avatar
Dion Haefner committed
672

673
          bcType = BCType::Neumann;
Dion Haefner's avatar
Dion Haefner committed
674 675 676 677
        }

        else if (BoundaryCondition::isDirichlet(bc))
        {
678
          bcType = BCType::Dirichlet;
Dion Haefner's avatar
Dion Haefner committed
679 680 681 682 683 684 685
        }

        else if (BoundaryCondition::isOther(bc))
          bcType = BCType::none;


        // update residual according to byType flag
686
        if (bcType == BCType::Neumann)
Dion Haefner's avatar
Dion Haefner committed
687 688
        {
          // flux is given parallel to gravity vector
689
          normal_flux *= std::abs( ig.centerUnitOuterNormal() * param->gravity() );
Dion Haefner's avatar
Dion Haefner committed
690 691

          // update residual (flux term)
692 693 694
          if constexpr (lopcase != LOPCase::RTVolume)
            for (unsigned int i = 0; i<lfsv_s.size(); i++)
              r_s.accumulate(lfsv_s,i, normal_flux * phiv_s[i] * factor);
Dion Haefner's avatar
Dion Haefner committed
695 696 697 698

          continue;
        }

699
        else if (bcType == BCType::Dirichlet)
Dion Haefner's avatar
Dion Haefner committed
700 701 702 703
        {
          // compute gradient of u
          Vector gradu_s(0.);
          for (unsigned int i = 0; i<lfsu_s.size(); i++)
704
            gradu_s.axpy(x_s(lfsu_s,i),gradphiu_s[i]);
Dion Haefner's avatar
Dion Haefner committed
705 706

          // update with gravity vector
707
          gradu_s -= param->gravity();
Dion Haefner's avatar
Dion Haefner committed
708 709 710 711 712

          // face normal vector
          const Domain normal = ig.unitOuterNormal(it.position());

          // jump relative to Dirichlet value
713
          const RF g = boundary->g(ig.intersection(),it.position(),time);
Dion Haefner's avatar
Dion Haefner committed
714 715
          const RF jump = u_s - g;

716 717 718 719 720
          // conductivity factors
          const RF relCond_s =
            conductivity_f_s(saturation_f_s(u_s)) / satCond_s;
          const RF relCond_n =
            conductivity_f_s(saturation_f_s(g)) / satCond_s;
721 722

          // penalty
Dion Haefner's avatar
Dion Haefner committed
723 724 725
          const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);

          // compute numerical flux
726
          RF numFlux = satCond_s*(penalty * jump - gradu_s * normal);
Dion Haefner's avatar
Dion Haefner committed
727

728 729
          // compute conductivity factor
          RF relCond;
730 731
          if (upwinding == RichardsUpwinding::semiUpwind
            || upwinding == RichardsUpwinding::fullUpwind)
732 733 734 735 736 737
          {
            if (numFlux > 0)
              relCond = relCond_s;
            else
              relCond = relCond_n;
          }
738
          else // (upwinding == RichardsUpwinding::none)
739 740 741 742
          {
            // harmonic average of conductivity factors
            relCond = 2 * relCond_s * relCond_n / (relCond_s + relCond_n);
          }
Dion Haefner's avatar
Dion Haefner committed
743 744


745 746 747 748 749 750 751 752 753
          if constexpr (lopcase != LOPCase::RTVolume)
          {
            // update residual (flux term)
            // diffusion term
            // consistency term
            // + penalty term
            for (unsigned int i = 0; i<lfsv_s.size(); i++)
              r_s.accumulate(lfsv_s,i, relCond * numFlux * phiv_s[i] * factor);
          }
Dion Haefner's avatar
Dion Haefner committed
754

755 756 757 758 759
          if constexpr (lopcase != LOPCase::RTSkeleton)
          {
            // update residual (symmetry term)
            // (non-)symmetric IP term
            // symmetry term
760 761 762 763

            // discrete gradient sign: Negative for lifting
            double dg_sign = (lopcase == LOPCase::RTVolume) ? -1. : 1.;

764 765 766
            for (unsigned int i = 0; i<lfsv_s.size(); i++)
              r_s.accumulate(lfsv_s,i, dg_sign * theta * (gradphiv_s[i] * normal) * relCond * satCond_s * jump * factor);
          }
Dion Haefner's avatar
Dion Haefner committed
767 768 769 770 771 772 773 774 775 776 777
          continue;
        }

        else if (bcType == BCType::none)
          continue;

        // The following implementation leads to some upwinding and nicer numerical behavior, but does not allow for water
        // stowage on the soil surface. Consider reimplementation if numerical problems occur with a limitedFlux boundary condition.

/*            if (BoundaryCondition::isLimitedFlux(bc))
        {
778 779 780
          const RF satCond_s   = param->K(*(ig.inside()),local_s);
          normal_flux = boundary->j(ig.intersection(),it.position(),time);
          normal_flux *= std::abs( ig.centerUnitOuterNormal() * param->gravity() );
Dion Haefner's avatar
Dion Haefner committed
781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820

          if (normal_flux < - satCond_s)
          {
            //std::cout << "- satCond " << - satCond_s << " flux " << normal_flux << std::endl;
            for (unsigned int i = 0; i < lfsv_s.size(); i++)
              r_s.accumulate(lfsv_s,i, - satCond_s * phi_s[i] * factor);
          }
          else if (normal_flux > satCond_s)
          {
            //std::cout << "  satCond " <<   satCond_s << " flux " << normal_flux << std::endl;
            for (unsigned int i = 0; i < lfsv_s.size(); i++)
              r_s.accumulate(lfsv_s,i,   satCond_s * phi_s[i] * factor);
          }
          else
          {
            for (unsigned int i = 0; i < lfsv_s.size(); i++)
              r_s.accumulate(lfsv_s,i, normal_flux * phi_s[i] * factor);
          }

          continue;
        }*/

      }
    }

  /**
   * @brief Volume integral depending only on test functions
   */
  template<typename EG, typename LFSV, typename R>
    void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
    {
      const int order = lfsv.finiteElement().localBasis().order();
      const int intorder = intorderadd + quadrature_factor * order;

      // get element geometry
      auto gt = eg.geometry();

      // loop over quadrature points
      for (const auto& it : quadratureRule(gt,intorder))
      {
821
        // evaluate position in element local coordinates
822
        const auto p_local = it.position();
Dion Haefner's avatar
Dion Haefner committed
823 824 825

        // evaluate values of shape functions (we assume Galerkin method lfsu = lfsv)
        std::vector<Scalar> phi(lfsv.size());
826
        lfsv.finiteElement().localBasis().evaluateFunction(p_local,phi);
Dion Haefner's avatar
Dion Haefner committed
827 828

        // scaled weight factor
829
        const RF factor = it.weight() * gt.integrationElement(p_local);
Dion Haefner's avatar
Dion Haefner committed
830

831
        const RF source = sourceTerm->q(eg.entity(),p_local,time) / gt.volume();
Dion Haefner's avatar
Dion Haefner committed
832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855

        // update residual
        for (unsigned int i = 0; i<lfsv.size(); i++)
          r.accumulate(lfsv,i, - source * phi[i] * factor);

        //std::cout << "lambda " << - source << std::endl;
      }
    }

  //! set operator time
  void setTime (RF t)
  {
    time = t;
  }

  //! set operator time and erase BC cache map
  void preStep (RF t, RF dt, int stages)
  {
    time = t;
  }
};

/**
 * @brief DG local operator for the temporal derivative of the Richards equation
856
 *
857
 * @ingroup LocalOperators
858
 * @ingroup RichardsModel
Dion Haefner's avatar
Dion Haefner committed
859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887
 */
template<typename Traits, typename Parameter, typename FEM, bool adjoint>
  class RichardsDGTemporalOperator
  : public Dune::PDELab::NumericalJacobianVolume<RichardsDGTemporalOperator<Traits, Parameter, FEM, adjoint> >,
  public Dune::PDELab::NumericalJacobianApplyVolume<RichardsDGTemporalOperator<Traits, Parameter, FEM, adjoint> >,
  public Dune::PDELab::FullSkeletonPattern,
  public Dune::PDELab::FullVolumePattern,
  public Dune::PDELab::LocalOperatorDefaultFlags,
  public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<typename Traits::RangeField>
{
  //typedef typename Parameter::Traits::GridTraits GridTraits;

  enum { dim = Traits::dim };

  typedef typename Traits::RangeField RF;
  typedef typename Traits::Scalar Scalar;
  typedef typename Traits::DomainField DF;
  typedef typename Traits::Domain Domain;

  public:

  // pattern assembly flags
  enum { doPatternVolume = true };

  // residual assembly flags
  enum { doAlphaVolume = true };

  private:

888
  std::shared_ptr<const Parameter> param; //!< class providing equation parameters
Dion Haefner's avatar
Dion Haefner committed
889 890 891 892 893 894 895 896 897
  RF time;
  const int intorderadd;
  const int quadrature_factor;

  public:

  /**
   * @brief Constructor
   */
898 899
  RichardsDGTemporalOperator (
    const Dune::ParameterTree& config,
900
    std::shared_ptr<const Parameter> param_,
901 902 903 904 905 906 907 908
    int intorderadd_ = 0,
    int quadrature_factor_ = 2
  ) :
    Dune::PDELab::NumericalJacobianVolume<RichardsDGTemporalOperator<Traits, Parameter, FEM, adjoint> >(1.e-7),
    param(param_),
    intorderadd(intorderadd_),
    quadrature_factor(quadrature_factor_)
  { }
Dion Haefner's avatar
Dion Haefner committed
909 910 911 912 913 914 915 916 917 918 919 920 921 922

  /**
   * @brief Volume integral depending on test and ansatz functions
   */
  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
    void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
    {
      // get polynomial degree
      const int order    = lfsu.finiteElement().localBasis().order();
      const int intorder = intorderadd + quadrature_factor * order;

      // get element geometry
      auto gt = eg.geometry();

923
      // bind parameterization and retrieve functions
924 925 926
      param->bind(eg.entity());
      const auto saturation_f = param->saturation_f();
      const auto water_content_f = param->water_content_f();
927

Dion Haefner's avatar
Dion Haefner committed
928 929 930
      // loop over quadrature points
      for (const auto& it : quadratureRule(gt,intorder))
      {
931
        // evaluate position in element local and global coordinates
932
        const auto p_local = it.position();
933

Dion Haefner's avatar
Dion Haefner committed
934 935
        // evaluate basis functions
        std::vector<Scalar> phi(lfsu.size());
936
        lfsu.finiteElement().localBasis().evaluateFunction(p_local, phi);
Dion Haefner's avatar
Dion Haefner committed
937 938 939 940

        // evaluate u
        RF u = 0.;
        for (unsigned int i = 0; i<lfsu.size(); i++)
941
          u += x(lfsu,i) * phi[i];
Dion Haefner's avatar
Dion Haefner committed
942

943
        // calculate water content from matric head
944
        const RF water_content = water_content_f(saturation_f(u));
Dion Haefner's avatar
Dion Haefner committed
945

946
        // integration factor
947
        const RF factor = it.weight() * gt.integrationElement(p_local);
Dion Haefner's avatar
Dion Haefner committed
948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969

        // update residual
        for (unsigned int i = 0; i<lfsv.size(); i++)
          r.accumulate(lfsv,i, water_content * phi[i] * factor);

      }
    }

  //! set operator time
  void setTime (RF t)
  {
    time = t;
  }


};

} // namespace Operator
} // namespace Dorie
} // namespace Dune


970
#endif // DUNE_DORIE_RICHARDS_OPERATOR_HH