local_operator_FV.hh 15 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
#ifndef DUNE_DORIE_RICHARDS_OPERATOR_FV_HH
#define DUNE_DORIE_RICHARDS_OPERATOR_FV_HH

#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>
13
#include <dune/dorie/model/richards/local_operator_cfg.hh>
14 15 16 17 18

namespace Dune {
namespace Dorie {
namespace Operator {

19 20 21 22 23 24 25 26 27 28 29 30 31 32
/*-------------------------------------------------------------------------*//**
 * @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
 */
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
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 };

55 56 57 58 59 60 61

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

72 73 74 75 76 77 78 79 80 81 82 83
   /*----------------------------------------------------------------------*//**
    * @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
84 85
    * @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
86 87 88 89 90 91 92
    *
    * @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
    */
93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
  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
118 119 120
    const auto geo_f = entity_f.geometry();
    const auto geo_i = entity_i.geometry();
    const auto geo_o = entity_o.geometry();
121 122 123 124 125

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

    // Get normal
126
    const auto normal_f = entity_f.centerUnitOuterNormal();
127

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

132 133 134 135 136
    // 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();

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

141 142 143 144 145 146
    // 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;

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

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

160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199
    // 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);
      if (dudn < 0.0) {
        // 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
      if (dudn < 0.0) {
        // inward flux
        cond_f = cond_o;
      }
      else {
        // outward flux
        cond_f = cond_i;
      }
    }
200 201

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

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

209 210 211 212 213 214 215 216
  /*-----------------------------------------------------------------------*//**
   * @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
217
   * @param[out] r_i     The view of the residual vector w.r.t lfsu_i
218 219 220 221 222 223 224
   *
   * @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
   */
225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245
  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
246
    const auto geo_f = entity_f.geometry();
247 248
    const auto& ref_el_f = referenceElement(geo_f);
    const auto& center_position_f = ref_el_f.position(0,0);
249 250 251 252 253 254 255 256

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

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

    // Retrieve boundary condition
257
    auto bc = _boundary->bc(entity_f,center_position_f,_time);
258 259

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

    if (BoundaryCondition::isDirichlet(bc))
    {
266 267
      const auto geo_i = entity_i.geometry();
      const auto center_position_i_g = geo_i.center();
268 269

      // Face center in global coordinates
270
      const auto center_position_f_g = geo_f.center();
271 272

      // Compute distance of these two points
273 274
      const auto position_diff = center_position_i_g - center_position_f_g;
      const auto distance = position_diff.two_norm();
275 276

      // Evaluate Dirichlet condition
277
      const auto g = _boundary->g(entity_f,center_position_f,_time);
278 279

      // Inside unknown value
280
      const RangeU u_i = x_i(lfsu_i,0);
281 282 283 284 285 286 287 288

      // 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
289
      dudn -= _param->gravity()*normal_f;
290 291

      // Water flux in normal direction w.r.t the intersection
292
      const auto water_flux_n = - cond_i*dudn;
293 294 295 296 297 298 299

      // 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
300
      RangeFieldU water_flux_f = _boundary->j(entity_f,center_position_f,_time);
301 302

      // Convert water flux into units of matric head flux
303
      water_flux_f *= std::abs( normal_f * _param->gravity() );
304
      // water_flux_f *= _param->gravity().two_norm();
305 306 307 308 309

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

310 311 312 313 314 315 316
  /*-----------------------------------------------------------------------*//**
   * @brief      Sets the time.
   *
   * @param[in]  t           time of type RangeField
   *
   * @tparam     RangeField  type of the range field
   */
317 318 319 320 321 322 323
  template<class RangeField>
  void setTime (RangeField t)
  {
    _time = t;
  }

private:
324
  /// Class providing equation parameters
325
  const std::shared_ptr<const Parameter> _param;
326
  /// Class providing boundary conditions
327
  const std::shared_ptr<const Boundary> _boundary;
328 329 330
  /// Upwinding setting
  const RichardsUpwinding _upwinding;
  /// Operator internal time
331 332 333
  double _time;
};

334 335 336 337 338 339 340 341 342 343 344
/**
 * @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
 */
345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365
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:

366 367 368 369 370
  /**
   * @brief      Constructor of RichardsFVTemporalOperator
   *
   * @param[in]  param     The parameter class
   */
371
  RichardsFVTemporalOperator(const std::shared_ptr<const Parameter>& param)
372 373 374 375
    : _param(param)
    , _time(0.)
  {}

376 377 378
  /*-----------------------------------------------------------------------*//**
   * @brief      Volume integral depending on test and ansatz functions
   *
379
   * @param[in]  eg    The entity of the grid
380 381 382
   * @param[in]  lfsu  The ansatz local function space
   * @param[in]  x     The coefficients of the lfsu
   * @param[in]  lfsv  The test local function space
383
   * @param[out] r     The view of the residual vector w.r.t lfsu
384 385 386 387 388 389 390
   *
   * @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
   */
391 392 393 394 395 396 397 398 399
  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;
400
    const RangeU u = x(lfsu,0);
401 402

    // entity geometry
403
    const auto geo = eg.geometry();
404 405

    // bind parameterization and retrieve functions
406 407 408
    _param->bind(eg.entity());
    const auto saturation_f = _param->saturation_f();
    const auto water_content_f = _param->water_content_f();
409 410 411 412 413 414 415 416

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

417
  /**
418 419
   * @brief      Sets the time.
   *
420
   * @param[in]  t     time of type RangeField
421
   *
422
   * @tparam     RF    Type of the range field
423 424 425 426 427 428 429 430
   */
  template<class RF>
  void setTime (RF t)
  {
    _time = t;
  }

private:
431
  /// class providing equation parameters
432
  const std::shared_ptr<const Parameter> _param;
433 434 435 436 437 438 439 440 441 442
  double _time;
};


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


#endif // DUNE_DORIE_RICHARDS_OPERATOR_FV_HH