local_operator_FV.hh 15.3 KB
Newer Older
1 2 3
#ifndef DUNE_DORIE_RICHARDS_OPERATOR_FV_HH
#define DUNE_DORIE_RICHARDS_OPERATOR_FV_HH

4 5
#include <dune/common/float_cmp.hh>

6 7 8 9 10 11 12 13 14
#include <dune/geometry/referenceelements.hh>

#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/dorie/common/typedefs.hh>
15
#include <dune/dorie/model/richards/local_operator_cfg.hh>
16 17 18 19 20

namespace Dune {
namespace Dorie {
namespace Operator {

21 22 23 24 25 26 27 28 29 30 31 32 33 34
/*-------------------------------------------------------------------------*//**
 * @brief      Spatial local operator for the Richards equation in a finite 
 *             volume scheme.
 * @details    It solves the spatial part of the Richards equation.
 * 
 * @author     Santiago Ospina De Los Ríos
 * @date       2019
 * @ingroup    LocalOperators
 * @ingroup    RichardsModel
 * 
 * @tparam     Parameter                 Type of the class providing parameters
 * @tparam     Boundary                  Type of the class providing boundary
 *                                       conditions
 */
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
template<class Parameter, class Boundary>
class RichardsFVSpatialOperator
  : public Dune::PDELab::NumericalJacobianVolume<
              RichardsFVSpatialOperator<Parameter,Boundary>>
  , public Dune::PDELab::NumericalJacobianSkeleton<
              RichardsFVSpatialOperator<Parameter,Boundary>>
  , public Dune::PDELab::NumericalJacobianBoundary<
              RichardsFVSpatialOperator<Parameter,Boundary>>
  , public Dune::PDELab::FullSkeletonPattern
  , public Dune::PDELab::FullVolumePattern
  , public Dune::PDELab::LocalOperatorDefaultFlags
  , public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>
{
public:
  // Pattern assembly flags
  enum { doPatternVolume = true };
  enum { doPatternSkeleton = true };

  // Residual assembly flags
  enum { doAlphaBoundary = true };
  enum { doAlphaSkeleton  = true };

57 58 59 60 61 62 63

  /**
   * @brief      Constructor of RichardsFVSpatialOperator
   *
   * @param[in]  param     The parameter class
   * @param[in]  boundary  The boundary class
   */
64
  RichardsFVSpatialOperator(
65
    const std::shared_ptr<const Parameter>& param,
66 67
    const std::shared_ptr<const Boundary>& boundary,
    const RichardsUpwinding upwinding = RichardsUpwinding::none
68 69
  ) : _param(param)
    , _boundary(boundary)
70
    , _upwinding(upwinding)
71 72 73
    , _time(0.)
  {}

74 75 76 77 78 79 80 81 82 83 84 85
   /*----------------------------------------------------------------------*//**
    * @brief      Skeleton integral depending on test and ansatz functions. Each
    *             face is only visited once since this method is symmetric
    *
    * @param[in]  ig      The intersection entity of the grid (inside + outside
    *                     entities)
    * @param[in]  lfsu_i  The inside ansatz local function space
    * @param[in]  x_i     The coefficients of the lfsu_i
    * @param[in]  lfsv_i  The inside test local function space
    * @param[in]  lfsu_o  The outside ansatz local function space
    * @param[in]  x_o     The coefficients of the lfsu_o
    * @param[in]  lfsv_o  The outside test local function space
86 87
    * @param[out] r_i     The view of the residual vector w.r.t lfsu_i
    * @param[out] r_o     The view of the residual vector w.r.t lfsu_o
88 89 90 91 92 93 94
    *
    * @tparam     IG      The type for ig
    * @tparam     LFSU    The type for lfsu_i and lfsu_o
    * @tparam     X       The type for x_i and x_o
    * @tparam     LFSV    The type for lfsv_i and lfsv_o
    * @tparam     R       The type for r_i and r_o
    */
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
  void alpha_skeleton(const IG& ig, 
    const LFSU& lfsu_i, const X& x_i, const LFSV& lfsv_i,
    const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o, 
    R& r_i, R& r_o) const
  {
    // Get local basis traits from local function space
    using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
                        Traits::LocalBasisType::Traits;

    // Get range, range field and gradient for trial space
    using RangeU = typename LocalBasisTraitsU::RangeType;

    // Assert that polynomial degree is always 0
    assert(lfsu_i.finiteElement().localBasis().order()==0);
    assert(lfsu_o.finiteElement().localBasis().order()==0);
    assert(lfsv_i.finiteElement().localBasis().order()==0);
    assert(lfsv_o.finiteElement().localBasis().order()==0);

    // Define entities for each frame of reference
    const auto& entity_f = ig.intersection();
    const auto& entity_i = ig.inside();
    const auto& entity_o = ig.outside();

    // Get geometries
120 121 122
    const auto geo_f = entity_f.geometry();
    const auto geo_i = entity_i.geometry();
    const auto geo_o = entity_o.geometry();
123 124 125 126 127

    // Get volume of entities
    const auto volume_f = geo_f.volume();

    // Get normal
128
    const auto normal_f = entity_f.centerUnitOuterNormal();
129

130
    // Entity centers in global coordinates
131 132
    const auto center_position_i_g = geo_i.center();
    const auto center_position_o_g = geo_o.center();
133

134 135 136 137 138
    // Distance between the two entity centers
    const auto center_position_diff = center_position_i_g
                                      - center_position_o_g;
    const auto distance = center_position_diff.two_norm();

139 140 141 142
    // Inside/outside solute value
    const RangeU u_i = x_i(lfsu_i,0);
    const RangeU u_o = x_o(lfsu_o,0);

143 144 145 146 147 148
    // Finite difference of u between the two entities
    RangeU dudn = (u_o-u_i) / distance;

    // Update gradient with gravity vector
    dudn -= _param->gravity() * normal_f;

149
    // bind parameterization and retrieve functions
150 151 152
    _param->bind(entity_i);
    const auto saturation_f_i = _param->saturation_f();
    const auto conductivity_f_i = _param->conductivity_f();
153
    
154 155 156
    _param->bind(entity_o);
    const auto saturation_f_o = _param->saturation_f();
    const auto conductivity_f_o = _param->conductivity_f();
157 158 159 160 161

    // Compute the conductivity
    const RangeU cond_i = conductivity_f_i(saturation_f_i(u_i));
    const RangeU cond_o = conductivity_f_o(saturation_f_o(u_o));

162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
    // Interface conductivity factor (Upwinding applies)
    RangeU cond_f(0.0);
    if (_upwinding == RichardsUpwinding::none)
    {
      // Harmonic average for interface conductivity
      cond_f = 2.0/(   1.0/(cond_i + 1E-30)
                     + 1.0/(cond_o + 1E-30) );
    }
    else if (_upwinding == RichardsUpwinding::semiUpwind)
    {
      // Harmonic average of saturated conductivities
      const RangeU cond_sat_i = conductivity_f_i(1.0);
      const RangeU cond_sat_o = conductivity_f_o(1.0);

      // Upwinding conductivity factor
      RangeU cond_upwind_factor(0.0);
178 179
      // NOTE: Sign seems inverted due to definition of 'water_flux_n' below.
      if (dudn > 0.0) {
180 181 182 183 184 185 186 187 188 189 190 191 192 193
        // inward flux
        cond_upwind_factor = cond_o / cond_sat_o;
      }
      else {
        // outward flux
        cond_upwind_factor = cond_i / cond_sat_i;
      }

      cond_f = 2.0 * cond_upwind_factor / (   1.0/(cond_sat_i + 1E-30)
                                            + 1.0/(cond_sat_o + 1E-30) );
    }
    else // RichardsFVUpwinding::fullUpwind
    {
      // No average: Use upwind conductivity
194
      if (dudn > 0.0) {
195 196 197 198 199 200 201 202
        // inward flux
        cond_f = cond_o;
      }
      else {
        // outward flux
        cond_f = cond_i;
      }
    }
203 204

    // Water flux in normal direction w.r.t the intersection
205
    const auto water_flux_n = - cond_f * dudn;
206 207 208 209 210 211

    // Symmetric contribution to residual on inside element
    r_i.accumulate(lfsv_i, 0,  water_flux_n*volume_f );
    r_o.accumulate(lfsv_o, 0, -water_flux_n*volume_f );
  }

212 213 214 215 216 217 218 219
  /*-----------------------------------------------------------------------*//**
   * @brief      Boundary integral depending on test and ansatz functions.
   *
   * @param[in]  ig      The intersection entity of the grid (inside + outside
   *                     entities)
   * @param[in]  lfsu_i  The inside ansatz local function space
   * @param[in]  x_i     The coefficients of the lfsu_i
   * @param[in]  lfsv_i  The inside test local function space
220
   * @param[out] r_i     The view of the residual vector w.r.t lfsu_i
221 222 223 224 225 226 227
   *
   * @tparam     IG      The type for ig
   * @tparam     LFSU    The type for lfsu_i
   * @tparam     X       The type for x_i
   * @tparam     LFSV    The type for lfsv_i
   * @tparam     R       The type for r_i
   */
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
  void alpha_boundary (const IG& ig, const LFSU& lfsu_i, const X& x_i,
                       const LFSV& lfsv_i, R& r_i) const
  {
    // Get local basis traits from local function space
    using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
                        Traits::LocalBasisType::Traits;

    // Get range, range field and gradient for trial space
    using RangeU = typename LocalBasisTraitsU::RangeType;
    using RangeFieldU = typename LocalBasisTraitsU::RangeFieldType;
    
    // Assert that polynomial degree is always 0
    assert(lfsu_i.finiteElement().localBasis().order()==0);
    assert(lfsv_i.finiteElement().localBasis().order()==0);

    // Get entities
    const auto& entity_f = ig.intersection(); 
    const auto& entity_i = ig.inside();

    // Get geometries
249
    const auto geo_f = entity_f.geometry();
250 251
    const auto& ref_el_f = referenceElement(geo_f);
    const auto& center_position_f = ref_el_f.position(0,0);
252 253 254 255 256 257 258 259

    // Face volume for integration
    const auto volume_f = geo_f.volume();

    // Normal vector of intersection
    const auto normal_f = ig.centerUnitOuterNormal();

    // Retrieve boundary condition
260
    auto bc = _boundary->bc(entity_f,center_position_f,_time);
261 262

    // bind parameterization and retrieve functions
263 264 265
    _param->bind(entity_i);
    const auto saturation_f_i = _param->saturation_f();
    const auto conductivity_f_i = _param->conductivity_f();
266

267
    if (BoundaryCondition::isDirichlet(bc) || BoundaryCondition::isOutflow(bc))
268
    {
269 270
      const auto geo_i = entity_i.geometry();
      const auto center_position_i_g = geo_i.center();
271 272

      // Face center in global coordinates
273
      const auto center_position_f_g = geo_f.center();
274 275

      // Compute distance of these two points
276 277
      const auto position_diff = center_position_i_g - center_position_f_g;
      const auto distance = position_diff.two_norm();
278 279

      // Evaluate Dirichlet condition
280
      const auto g = _boundary->g(entity_f,center_position_f,_time);
281 282

      // Inside unknown value
283
      const RangeU u_i = x_i(lfsu_i,0);
284 285 286 287 288 289 290 291

      // Compute the conductivity
      const RangeU cond_i = conductivity_f_i(saturation_f_i(u_i));

      // Finite difference of u between the two entities
      RangeU dudn = (g-u_i)/distance;

      // Update gradient with gravity vector
292
      dudn -= _param->gravity()*normal_f;
293 294

      // Water flux in normal direction w.r.t the intersection
295
      const auto water_flux_n = - cond_i*dudn;
296

297 298 299 300 301
      // No contribution for inward flux
      if (BoundaryCondition::isOutflow(bc)
          and Dune::FloatCmp::lt(double(water_flux_n), 0.0))
        return;

302 303 304 305 306 307
      // Contribution to residual from Dirichlet boundary
      r_i.accumulate(lfsv_i, 0, water_flux_n*volume_f);
    }
    else if (BoundaryCondition::isNeumann(bc))
    {
      // Contribution to residual from Neumann boundary
308
      RangeFieldU water_flux_f = _boundary->j(entity_f,center_position_f,_time);
309 310

      // Convert water flux into units of matric head flux
311
      water_flux_f *= std::abs( normal_f * _param->gravity() );
312
      // water_flux_f *= _param->gravity().two_norm();
313 314 315 316 317

      r_i.accumulate(lfsv_i, 0, water_flux_f*volume_f );
    }
  }

318 319 320 321 322 323 324
  /*-----------------------------------------------------------------------*//**
   * @brief      Sets the time.
   *
   * @param[in]  t           time of type RangeField
   *
   * @tparam     RangeField  type of the range field
   */
325 326 327 328 329 330 331
  template<class RangeField>
  void setTime (RangeField t)
  {
    _time = t;
  }

private:
332
  /// Class providing equation parameters
333
  const std::shared_ptr<const Parameter> _param;
334
  /// Class providing boundary conditions
335
  const std::shared_ptr<const Boundary> _boundary;
336 337 338
  /// Upwinding setting
  const RichardsUpwinding _upwinding;
  /// Operator internal time
339 340 341
  double _time;
};

342 343 344 345 346 347 348 349 350 351 352
/**
 * @brief      Temporal local operator Richards equation in a finite volume
 *             scheme.
 * @details    It solves the temporal part of the Richards equation.
 * @author     Santiago Ospina De Los Ríos
 * @date       2018
 * @ingroup    LocalOperators
 * @ingroup    RichardsModel
 *
 * @tparam     Parameter                 Type of the class providing parameters
 */
353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373
template<class Parameter>
  class RichardsFVTemporalOperator
  : public Dune::PDELab::NumericalJacobianVolume<
              RichardsFVTemporalOperator<Parameter>>
  , public Dune::PDELab::NumericalJacobianApplyVolume<
              RichardsFVTemporalOperator<Parameter>>
  , public Dune::PDELab::FullSkeletonPattern
  , public Dune::PDELab::FullVolumePattern
  , public Dune::PDELab::LocalOperatorDefaultFlags
  , public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>
{
public:

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

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

public:

374 375 376 377 378
  /**
   * @brief      Constructor of RichardsFVTemporalOperator
   *
   * @param[in]  param     The parameter class
   */
379
  RichardsFVTemporalOperator(const std::shared_ptr<const Parameter>& param)
380 381 382 383
    : _param(param)
    , _time(0.)
  {}

384 385 386
  /*-----------------------------------------------------------------------*//**
   * @brief      Volume integral depending on test and ansatz functions
   *
387
   * @param[in]  eg    The entity of the grid
388 389 390
   * @param[in]  lfsu  The ansatz local function space
   * @param[in]  x     The coefficients of the lfsu
   * @param[in]  lfsv  The test local function space
391
   * @param[out] r     The view of the residual vector w.r.t lfsu
392 393 394 395 396 397 398
   *
   * @tparam     EG    The type for eg
   * @tparam     LFSU  The type for lfsu
   * @tparam     X     The type for x
   * @tparam     LFSV  The type for lfsv
   * @tparam     R     The type for r
   */
399 400 401 402 403 404 405 406 407
  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 local basis traits from local function space
    using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
                        Traits::LocalBasisType::Traits;
    // Get range for trial space
    using RangeU = typename LocalBasisTraitsU::RangeType;
408
    const RangeU u = x(lfsu,0);
409 410

    // entity geometry
411
    const auto geo = eg.geometry();
412 413

    // bind parameterization and retrieve functions
414 415 416
    _param->bind(eg.entity());
    const auto saturation_f = _param->saturation_f();
    const auto water_content_f = _param->water_content_f();
417 418 419 420 421 422 423 424

    // Calculate water content from matric head
    const RangeU water_content = water_content_f(saturation_f(u));

    // update residual
    r.accumulate(lfsv ,0 , water_content*geo.volume());
  }

425
  /**
426 427
   * @brief      Sets the time.
   *
428
   * @param[in]  t     time of type RangeField
429
   *
430
   * @tparam     RF    Type of the range field
431 432 433 434 435 436 437 438
   */
  template<class RF>
  void setTime (RF t)
  {
    _time = t;
  }

private:
439
  /// class providing equation parameters
440
  const std::shared_ptr<const Parameter> _param;
441 442 443 444 445 446 447 448 449 450
  double _time;
};


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


#endif // DUNE_DORIE_RICHARDS_OPERATOR_FV_HH