local_operator.hh 31.6 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 6 7 8 9 10 11 12
#include <string>
#include <tuple>
#include <map>
#include <vector>

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

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

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

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

24
#include <dune/dorie/common/typedefs.hh>
Dion Haefner's avatar
Dion Haefner committed
25 26 27 28 29 30 31 32 33

namespace Dune {
namespace Dorie {
namespace Operator {

//! Interior penalty type
struct RichardsDGMethod
{
  enum Type { NIPG, //!< Non-symmetric interior penalty
34 35 36
              SIPG, //!< Symmetric interior penalty (default)
              OOB, //!< Oden, Babuska, Baumann (no penalty term)
              IIP //!< Incomplete interior penalty (no symmetry)
Dion Haefner's avatar
Dion Haefner committed
37 38 39
            };
};

40 41
//! Upwinding type
struct RichardsDGUpwinding
Dion Haefner's avatar
Dion Haefner committed
42
{
43 44 45 46 47
  enum Type {
    fullUpwind, //!< Upwind flux and conductiviy
    semiUpwind, //!< Only upwind conductivity factor
    none //!< No upwinding
  };
Dion Haefner's avatar
Dion Haefner committed
48 49
};

50 51 52 53 54 55 56 57 58
//! Gradient weighting
struct RichardsDGWeights
{
  enum Type {
    weightsOn, //!< harmonic weighting of flux contributions
    weightsOff //!< no weighting of flux contributions
  };
};

Dion Haefner's avatar
Dion Haefner committed
59 60 61
//! Operator internal BC handling
struct BCType
{
62
  enum Type { Neumann, //!< Fixed flux at boundary
63
              Dirichlet, //!< Fixed matric head at boundary
Dion Haefner's avatar
Dion Haefner committed
64 65 66 67
              none //!< No BC specified
            };
};

68 69 70 71 72 73 74 75 76 77
/// Read experimental operator settings from a parameter file
/* \return Tuple of RichardsDGMethod, RichardsDGUpwinding, and
 *    RichardsDGWeights to be inserted into the local operator constructor.
 */
inline auto read_experimental_operator_settings(
  const Dune::ParameterTree &inifile)
    -> std::tuple<Operator::RichardsDGMethod::Type,
                  Operator::RichardsDGUpwinding::Type,
                  Operator::RichardsDGWeights::Type>
{
78 79 80
  const auto log = Dorie::get_logger(log_richards);
  log->debug("Reading experimental DG operator settings:");

81 82 83 84 85
  RichardsDGMethod::Type method;
  RichardsDGUpwinding::Type upwinding;
  RichardsDGWeights::Type weights;

  const auto method_str = inifile.get<std::string>("dg.experimental.method");
86
  log->debug("  DG method: {}", method_str);
87 88 89 90 91 92 93 94
  if (method_str == "SIPG")
    method = RichardsDGMethod::SIPG;
  else if (method_str == "NIPG")
    method = RichardsDGMethod::NIPG;
  else if (method_str == "OOB")
    method = RichardsDGMethod::OOB;
  else if (method_str == "IIP")
    method = RichardsDGMethod::IIP;
95 96 97 98
  else {
    log->error("Unknown DG method: {}", method_str);
    DUNE_THROW(Dune::IOError, "Unknown DG method");
  }
99 100 101

  const auto upwinding_str
    = inifile.get<std::string>("dg.experimental.upwinding");
102
  log->debug("  Upwinding scheme: {}", upwinding_str);
103 104 105 106 107 108
  if (upwinding_str == "none")
    upwinding = RichardsDGUpwinding::none;
  else if (upwinding_str == "semiUpwind")
    upwinding = RichardsDGUpwinding::semiUpwind;
  else if (upwinding_str == "fullUpwind")
    upwinding = RichardsDGUpwinding::fullUpwind;
109 110 111 112
  else {
    log->error("Unknown upwinding scheme: {}", upwinding_str);
    DUNE_THROW(Dune::IOError, "Unknown upwinding scheme");
  }
113 114

  const auto weights_str = inifile.get<bool>("dg.experimental.weights");
115
  log->debug("  Harmonic weights at interfaces: {}", weights_str);
116 117 118 119 120 121 122 123
  if (weights_str)
    weights = RichardsDGWeights::weightsOn;
  else
    weights = RichardsDGWeights::weightsOff;

  return std::make_tuple(method, upwinding, weights);
}

Dion Haefner's avatar
Dion Haefner committed
124 125 126 127 128 129 130 131 132 133 134
/// 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.)
135
 *
136
 * @ingroup LocalOperators
Dion Haefner's avatar
Dion Haefner committed
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 164 165 166
 */
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 };
  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;

  const Parameter&  param; //!< class providing equation parameters
  const Boundary&   boundary; //!< class providing boundary conditions
  const SourceTerm& sourceTerm; //!< class providing source term information
167 168 169
  const RichardsDGMethod::Type method; //!< Interior penalty type
  const RichardsDGUpwinding::Type upwinding; //!< DG method weights switch
  const RichardsDGWeights::Type weights; //!< Gradient weighting
Dion Haefner's avatar
Dion Haefner committed
170 171 172

  const int intorderadd; //!< integration order addend
  const int quadrature_factor; //!< factor to FEM integration order
173 174 175
  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
176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195

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

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
196
   *  \param method_ DG method type \see RichardsDGMethod::Type
197
   *  \param upwinding_ Upwinding type \see RichardsDGUpwinding::Type
Dion Haefner's avatar
Dion Haefner committed
198 199 200 201
   *  \param weights_ DG weigths switch \see RichardsDGWeights::Type
   *  \param intorderadd_ Addend to integration order
   *  \param quadrature_factor_ Factor to FEM integration order
   */
202 203 204 205 206 207 208 209 210 211 212 213
  RichardsDGSpatialOperator (
    const Dune::ParameterTree& config,
    const Parameter& param_,
    const Boundary& boundary_,
    const SourceTerm& sourceTerm_,
    const RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
    const RichardsDGUpwinding::Type upwinding_ = RichardsDGUpwinding::none,
    const RichardsDGWeights::Type weights_ = RichardsDGWeights::weightsOn,
    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
214 215
    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),
216 217 218 219 220 221 222 223
    param(param_),
    boundary(boundary_),
    sourceTerm(sourceTerm_),
    method(method_),
    upwinding(upwinding_),
    weights(weights_),
    intorderadd(intorderadd_),
    quadrature_factor(quadrature_factor_),
224
    penalty_factor(config.get<RF>("numerics.penaltyFactor")),
225
    time(0.0),
Dion Haefner's avatar
Dion Haefner committed
226
    cache(20)
Dion Haefner's avatar
Dion Haefner committed
227
  {
228 229 230 231 232
    theta = 0.; // IIP
    if (method == RichardsDGMethod::NIPG || method == RichardsDGMethod::OOB)
      theta = 1.;
    else if (method == RichardsDGMethod::SIPG)
      theta = -1.;
233

234 235
    if (method == RichardsDGMethod::OOB)
      penalty_factor = 0.;
Dion Haefner's avatar
Dion Haefner committed
236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257
  }

  /// 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
    {
      const int order = lfsu.finiteElement().localBasis().order();
      const int intorder = intorderadd + quadrature_factor * order;

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

258 259 260 261 262
      // bind parameterization and retrieve functions
      param.bind(eg.entity());
      const auto saturation_f = param.saturation_f();
      const auto conductivity_f = param.conductivity_f();

Dion Haefner's avatar
Dion Haefner committed
263 264 265
      // loop over quadrature points
      for (const auto& it : quadratureRule(gt,intorder))
      {
266
        // evaluate position in element local and global coordinates
267
        const auto p_local = it.position();
268

Dion Haefner's avatar
Dion Haefner committed
269
        // evaluate basis functions
270
        auto& phi = cache[order].evaluateFunction(p_local,lfsu.finiteElement().localBasis());
Dion Haefner's avatar
Dion Haefner committed
271 272 273 274 275 276 277

        // evaluate u
        RF u = 0.;
        for (unsigned int i = 0; i<lfsu.size(); i++)
          u += x(lfsu,i) * phi[i];

        // evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
278
        const auto& js = cache[order].evaluateJacobian(p_local,lfsu.finiteElement().localBasis());
Dion Haefner's avatar
Dion Haefner committed
279 280

        // transform gradients to real element
281
        const auto jac = gt.jacobianInverseTransposed(p_local);
Dion Haefner's avatar
Dion Haefner committed
282 283 284 285 286 287 288 289 290 291 292
        std::vector<Vector> gradphi(lfsu.size());
        for (unsigned int i = 0; i<lfsu.size(); i++)
          jac.mv(js[i][0],gradphi[i]);

        // compute gradient of u
        Vector gradu(0.);
        for (unsigned int i = 0; i<lfsu.size(); i++)
          gradu.axpy(x(lfsu,i),gradphi[i]);

        gradu -= param.gravity();

293 294
        // compute conductivity
        const RF cond = conductivity_f(saturation_f(u));
Dion Haefner's avatar
Dion Haefner committed
295 296

        // integration factor
297
        const RF factor = it.weight() * gt.integrationElement(p_local);
Dion Haefner's avatar
Dion Haefner committed
298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341

        // update residual
        for (unsigned int i = 0; i<lfsv.size(); i++)
        {
          r.accumulate(lfsv,i, cond * (gradu*gradphi[i]) * factor);
        }
      }
    }

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

342 343 344 345 346 347
      // 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();

      // cache saturated conductivity
348
      const RF satCond_s = conductivity_f_s(1.0);
349 350 351 352 353

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

354
      const RF satCond_n = conductivity_f_n(1.0);
355

Dion Haefner's avatar
Dion Haefner committed
356 357 358 359 360 361 362
      // 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
363 364
        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
365

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

Dion Haefner's avatar
Dion Haefner committed
369
        // evaluate basis functions
370 371
        const auto& phi_s = cache[order_s].evaluateFunction(p_local_s,lfsu_s.finiteElement().localBasis());
        const auto& phi_n = cache[order_n].evaluateFunction(p_local_n,lfsu_n.finiteElement().localBasis());
Dion Haefner's avatar
Dion Haefner committed
372 373 374 375 376 377 378 379 380

        // 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)
381 382
        const auto& gradphi_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
        const auto& gradphi_n = cache[order_n].evaluateJacobian(p_local_n,lfsu_n.finiteElement().localBasis());
Dion Haefner's avatar
Dion Haefner committed
383 384 385 386 387

        // transform gradients to real element
        Tensor jac;
        std::vector<Vector> tgradphi_s(lfsu_s.size());
        std::vector<Vector> tgradphi_n(lfsu_n.size());
388
        jac = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
Dion Haefner's avatar
Dion Haefner committed
389 390
        for (unsigned int i = 0; i<lfsu_s.size(); i++)
          jac.mv(gradphi_s[i][0],tgradphi_s[i]);
391
        jac = ig.outside().geometry().jacobianInverseTransposed(p_local_n);
Dion Haefner's avatar
Dion Haefner committed
392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408
        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]);

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

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

409 410 411
        // 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));
412

Dion Haefner's avatar
Dion Haefner committed
413 414
        // compute weights
        RF omega_s = 0.5, omega_n = 0.5;
415
        if (weights == RichardsDGWeights::weightsOn)
416
        {
417 418
          if (upwinding == RichardsDGUpwinding::none)
          {
419 420
            omega_s = cond_n / ((cond_s) + (cond_n));
            omega_n = cond_s / ((cond_s) + (cond_n));
421 422 423 424 425 426 427
          }
          else if (upwinding == RichardsDGUpwinding::semiUpwind
            || upwinding == RichardsDGUpwinding::fullUpwind)
          {
            omega_s = satCond_n / (satCond_s + satCond_n);
            omega_n = satCond_s / (satCond_s + satCond_n);
          }
Dion Haefner's avatar
Dion Haefner committed
428 429 430 431 432
        }

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

433
        // penalty factor
Dion Haefner's avatar
Dion Haefner committed
434 435
        const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);

436 437 438 439
        // 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
440
        // compute numerical flux
441 442
        RF numFlux_s = satCond_s * ( - (gradu_s * normal) + penalty * jump);
        RF numFlux_n = satCond_n * ( - (gradu_n * normal) + penalty * jump);
443
        if (upwinding == RichardsDGUpwinding::none)
444 445 446 447 448
        {
          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
449

450
        // compute upwind quantities
451
        RF relCond_upwind=0.0, satCond_upwind=0.0, numFlux_upwind=0.0;
452 453 454
        if (upwinding == RichardsDGUpwinding::semiUpwind
          || upwinding == RichardsDGUpwinding::fullUpwind)
        {
455
          RF cond_upwind;
456 457
          if (numFlux > 0)
          {
458
            cond_upwind = conductivity_f_s(saturation_f_s(u_s));
459 460
            satCond_upwind = satCond_s;
            numFlux_upwind = numFlux_s;
461 462 463
          }
          else
          {
464
            cond_upwind = conductivity_f_n(saturation_f_n(u_n));
465
            satCond_upwind = satCond_n;
466
            numFlux_upwind = numFlux_n;
467
          }
468 469

          relCond_upwind = cond_upwind / satCond_upwind;
470
        }
Dion Haefner's avatar
Dion Haefner committed
471 472 473 474 475

        // update residual (flux term)
        // diffusion term
        // consistency term
        // + penalty term
476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492
        if (upwinding == RichardsDGUpwinding::none)
        {
          for (unsigned int i = 0; i < lfsv_s.size(); i++)
            r_s.accumulate(lfsv_s, i,   numFlux * phi_s[i] * factor);
          for (unsigned int i = 0; i < lfsv_n.size(); i++)
            r_n.accumulate(lfsv_n, i, - numFlux * phi_n[i] * factor);
        }
        else if (upwinding == RichardsDGUpwinding::semiUpwind)
        {
          for (unsigned int i = 0; i<lfsv_s.size(); i++)
            r_s.accumulate(lfsv_s,i,   relCond_upwind * numFlux * phi_s[i] * factor);
          for (unsigned int i = 0; i<lfsv_n.size(); i++)
            r_n.accumulate(lfsv_n,i, - relCond_upwind * numFlux * phi_n[i] * factor);
        }
        else if (upwinding == RichardsDGUpwinding::fullUpwind)
        {
          for (unsigned int i = 0; i < lfsv_s.size(); i++)
493
            r_s.accumulate(lfsv_s, i,   relCond_upwind * numFlux_upwind * phi_s[i] * factor);
494
          for (unsigned int i = 0; i < lfsv_n.size(); i++)
495
            r_n.accumulate(lfsv_n, i, - relCond_upwind * numFlux_upwind * phi_n[i] * factor);
496
        }
Dion Haefner's avatar
Dion Haefner committed
497 498 499 500

        // update residual (symmetry term)
        // (non-)symmetric IP term
        // symmetry term
501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521
        if (upwinding == RichardsDGUpwinding::none)
        {
          for (unsigned int i = 0; i < lfsv_s.size(); i++)
            r_s.accumulate(lfsv_s, i, theta * omega_s * (tgradphi_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, theta * omega_n * (tgradphi_n[i] * normal) * relCond_n * satCond_n * jump * factor);
        }
        else if (upwinding == RichardsDGUpwinding::semiUpwind)
        {
          for (unsigned int i = 0; i<lfsv_s.size(); i++)
            r_s.accumulate(lfsv_s,i, theta * omega_s * (tgradphi_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, theta * omega_n * (tgradphi_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
        }
        else if (upwinding == RichardsDGUpwinding::fullUpwind)
        {
          for (unsigned int i = 0; i < lfsv_s.size(); i++)
            r_s.accumulate(lfsv_s, i, theta * (tgradphi_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, theta * (tgradphi_n[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
        }
Dion Haefner's avatar
Dion Haefner committed
522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544
      }
    }

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

545 546 547 548 549 550
      // 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();

      // cache saturated conductivity
551
      const RF satCond_s = conductivity_f_s(1.0);
552

Dion Haefner's avatar
Dion Haefner committed
553 554 555 556
      // loop over quadrature points and integrate normal flux
      for (const auto& it : quadratureRule(gtface,intorder))
      {
        // position of quadrature point in local coordinates of elements
557
        const auto p_local_s = ig.geometryInInside().global(it.position());
Dion Haefner's avatar
Dion Haefner committed
558 559

        // evaluate basis functions
560
        const auto& phi_s = cache[order_s].evaluateFunction(p_local_s,lfsu_s.finiteElement().localBasis());
Dion Haefner's avatar
Dion Haefner committed
561 562 563 564 565 566 567 568 569

        // 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++)
          u_s += x_s(lfsu_s,i) * phi_s[i];

570
        // query the type of boundary condition
Dion Haefner's avatar
Dion Haefner committed
571
        // set the internal bcType flag
572
        const auto bc = boundary.bc(ig.intersection(),it.position(),time);
Dion Haefner's avatar
Dion Haefner committed
573 574 575 576 577 578 579 580
        auto bcType = BCType::none;
        RF normal_flux;

        if (BoundaryCondition::isNeumann(bc))
        {
          // evaluate flux boundary condition
          normal_flux = boundary.j(ig.intersection(),it.position(),time);

581
          bcType = BCType::Neumann;
Dion Haefner's avatar
Dion Haefner committed
582 583 584 585
        }

        else if (BoundaryCondition::isDirichlet(bc))
        {
586
          bcType = BCType::Dirichlet;
Dion Haefner's avatar
Dion Haefner committed
587 588 589 590 591 592 593
        }

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


        // update residual according to byType flag
594
        if (bcType == BCType::Neumann)
Dion Haefner's avatar
Dion Haefner committed
595 596 597 598 599 600 601 602 603 604 605
        {
          // flux is given parallel to gravity vector
          normal_flux *= std::abs( ig.centerUnitOuterNormal() * param.gravity() );

          // update residual (flux term)
          for (unsigned int i = 0; i<lfsv_s.size(); i++)
            r_s.accumulate(lfsv_s,i, normal_flux * phi_s[i] * factor);

          continue;
        }

606
        else if (bcType == BCType::Dirichlet)
Dion Haefner's avatar
Dion Haefner committed
607 608
        {
          // evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
609
          const auto& gradphi_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
Dion Haefner's avatar
Dion Haefner committed
610 611

          // transform gradients to real element
612
          const auto jac = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
Dion Haefner's avatar
Dion Haefner committed
613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631
          std::vector<Vector> tgradphi_s(lfsu_s.size());
          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]);

          // update with gravity vector
          gradu_s -= param.gravity();

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

          // jump relative to Dirichlet value
          const RF g = boundary.g(ig.intersection(),it.position(),time);
          const RF jump = u_s - g;

632 633 634 635 636
          // 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;
637 638

          // penalty
Dion Haefner's avatar
Dion Haefner committed
639 640 641 642 643
          const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);

          // compute numerical flux
          const RF numFlux = satCond_s * ( - (gradu_s * normal) + penalty * jump);

644 645 646 647 648 649 650 651 652 653
          // compute conductivity factor
          RF relCond;
          if (upwinding == RichardsDGUpwinding::semiUpwind
            || upwinding == RichardsDGUpwinding::fullUpwind)
          {
            if (numFlux > 0)
              relCond = relCond_s;
            else
              relCond = relCond_n;
          }
654
          else // (upwinding == RichardsDGUpwinding::none)
655 656 657 658
          {
            // harmonic average of conductivity factors
            relCond = 2 * relCond_s * relCond_n / (relCond_s + relCond_n);
          }
Dion Haefner's avatar
Dion Haefner committed
659 660 661 662 663 664

          // update residual (flux term)
          // diffusion term
          // consistency term
          // + penalty term
          for (unsigned int i = 0; i<lfsv_s.size(); i++)
665
            r_s.accumulate(lfsv_s,i, relCond * numFlux * phi_s[i] * factor);
Dion Haefner's avatar
Dion Haefner committed
666 667 668 669 670

          // update residual (symmetry term)
          // (non-)symmetric IP term
          // symmetry term
          for (unsigned int i = 0; i<lfsv_s.size(); i++)
671
            r_s.accumulate(lfsv_s,i, theta * (tgradphi_s[i] * normal) * relCond * satCond_s * jump * factor);
Dion Haefner's avatar
Dion Haefner committed
672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726

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

          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))
      {
727
        // evaluate position in element local coordinates
728
        const auto p_local = it.position();
Dion Haefner's avatar
Dion Haefner committed
729 730 731

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

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

737
        const RF source = sourceTerm.q(eg.entity(),p_local,time) / gt.volume();
Dion Haefner's avatar
Dion Haefner committed
738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761

        // 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
762
 *
763
 * @ingroup LocalOperators
Dion Haefner's avatar
Dion Haefner committed
764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802
 */
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:

  const Parameter& param; //!< class providing equation parameters
  RF time;
  const int intorderadd;
  const int quadrature_factor;

  public:

  /**
   * @brief Constructor
   */
803 804 805 806 807 808 809 810 811 812 813
  RichardsDGTemporalOperator (
    const Dune::ParameterTree& config,
    const Parameter& param_,
    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
814 815 816 817 818 819 820 821 822 823 824 825 826 827

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

828 829 830 831 832
      // bind parameterization and retrieve functions
      param.bind(eg.entity());
      const auto saturation_f = param.saturation_f();
      const auto water_content_f = param.water_content_f();

Dion Haefner's avatar
Dion Haefner committed
833 834 835
      // loop over quadrature points
      for (const auto& it : quadratureRule(gt,intorder))
      {
836
        // evaluate position in element local and global coordinates
837
        const auto p_local = it.position();
838

Dion Haefner's avatar
Dion Haefner committed
839 840
        // evaluate basis functions
        std::vector<Scalar> phi(lfsu.size());
841
        lfsu.finiteElement().localBasis().evaluateFunction(p_local, phi);
Dion Haefner's avatar
Dion Haefner committed
842 843 844 845

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

848
        // calculate water content from matric head
849
        const RF water_content = water_content_f(saturation_f(u));
Dion Haefner's avatar
Dion Haefner committed
850

851
        // integration factor
852
        const RF factor = it.weight() * gt.integrationElement(p_local);
Dion Haefner's avatar
Dion Haefner committed
853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874

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


875
#endif // DUNE_DORIE_RICHARDS_OPERATOR_HH