The TS-GitLab will have to shut down towards the end of the year — please think about migrating your projects to GitLab.com or GitHub.
(This is still a very early message, meant to keep you informed. There will be more urgent ones in the future.)

local_operator_FV.hh 21 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_DORIE_TRANSPORT_OPERATOR_HH
#define DUNE_DORIE_TRANSPORT_OPERATOR_HH

#include <vector>

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

15
#include <dune/dorie/common/typedefs.hh>
16 17 18 19 20 21 22 23

namespace Dune {
namespace Dorie {
namespace Operator {

/*-------------------------------------------------------------------------*//**
 * @brief      Spatial local operator for the transport equation in unsaturated
 *             media in a finite volume scheme.
24
 * @details    It solves the spatial part of the transport equation:
Lukas Riedel's avatar
Lukas Riedel committed
25 26
 *
 *             @f{eqnarray*}{
27 28
 *             \partial_t[\theta C_w] +
 *             \nabla\cdot [\textbf{j}_w C_w] -
Lukas Riedel's avatar
Lukas Riedel committed
29
 *             \nabla [\theta \mathsf{D}_{\text{eff}}\nabla C_w]&=0 &\qquad \text{in }
30
 *             \Omega\\
Lukas Riedel's avatar
Lukas Riedel committed
31
 *             C_w &= g &\qquad \text{on } \Gamma_D
32
 *             \subseteq\partial\Omega\\
Lukas Riedel's avatar
Lukas Riedel committed
33
 *             \nabla C_w \cdot \textbf{n} &= \textbf{j}_{\scriptscriptstyle
34
 *             C_w}& \qquad \text{on } \Gamma_N =\partial\Omega \backslash
Lukas Riedel's avatar
Lukas Riedel committed
35 36 37
 *             \Gamma_D
 *             @f}
 *
38 39 40
 * @author     Santiago Ospina De Los Ríos
 * @date       2018
 * @ingroup    LocalOperators
41
 * @ingroup    TransportModel
42
 *
43
 * @todo Use diffusion coefficient from parameterization
44 45 46 47 48 49
 *
 * @tparam     Boundary        Type of the class providing boundary conditions
 * @tparam     GFWaterFlux     Type of a grid function which provides the water
 *                             flux
 * @tparam     GFWaterContent  Type of a grid function which provides the water
 *                             content
50
 */
51
template<class Boundary, class GFWaterFlux, class GFWaterContent>
52 53 54 55
class TransportFVSpatialOperator
    : public Dune::PDELab::NumericalJacobianVolume<
                        TransportFVSpatialOperator<
                            Boundary,
56 57
                            GFWaterFlux,
                            GFWaterContent
58
                        >
Santiago Ospina's avatar
Santiago Ospina committed
59
                    > 
60 61 62
    , public Dune::PDELab::NumericalJacobianSkeleton<
                        TransportFVSpatialOperator<
                            Boundary,
63 64
                            GFWaterFlux,
                            GFWaterContent
65 66 67 68 69
                        >
                    >
    , public Dune::PDELab::NumericalJacobianBoundary<
                        TransportFVSpatialOperator<
                            Boundary,
70 71
                            GFWaterFlux,
                            GFWaterContent
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
                        >
                    >
    , 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 };
87

88
private:
89
    static_assert(std::is_same<
90 91 92 93
            typename GFWaterFlux::Traits::GridViewType,
            typename GFWaterContent::Traits::GridViewType>::value, 
        "TransportFVSpatialOperator: GFWaterFlux and"
        "GFWaterContent has to use the same grid view.");
94

95 96 97
    // Define range type for water flux and water content
    using WaterFlux = typename GFWaterFlux::Traits::RangeType;
    using WaterContent = typename GFWaterContent::Traits::RangeType;
98 99 100

public:

101 102 103
    /**
     * @brief      Constructor of TransportFVSpatialOperator
     *
104 105
     * @param      boundary          The boundary terms
     * @param[in]  diff_coeff        The diffusion coefficient
106
     */
107
    TransportFVSpatialOperator(
108
        std::shared_ptr<Boundary> boundary, 
109
        double diff_coeff
110
    )   : _boundary(boundary)
111
        , _time(0.)
112
        , _diff_coeff(diff_coeff)
113
    { }
114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141

    /*---------------------------------------------------------------------*//**
     * @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
     * @param      r_i     The view of the residual vector w.r.t lfsu_i
     * @param      r_o     The view of the residual vector w.r.t lfsu_o
     *
     * @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
     */
    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
    {
142 143 144
        // Get entity dim
        // const int dim = IG::Entity::dimension;

145 146
        // Get local basis traits from local function space
        using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
147 148 149
                                         Traits::LocalBasisType::Traits;

        // Get range, range field and gradient for trial space
150
        using RangeU = typename LocalBasisTraitsU::RangeType;
151 152
        // using RangeFieldU = typename LocalBasisTraitsU::RangeFieldType;
        // using GradientU = Dune::FieldVector<RangeFieldU,dim>;
153 154 155 156 157 158 159 160 161

        // 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();
162 163
        const auto& entity_i = ig.inside();
        const auto& entity_o = ig.outside();
164 165

        // Get geometries
166
        auto geo_f = entity_f.geometry();
167 168
        auto geo_i = entity_i.geometry();
        auto geo_o = entity_o.geometry();
169

170
        // Get geometry of intersection in local coord. of neighbor entities
171
        auto geo_in_i = entity_f.geometryInInside();
172

173
        // Get volume of entity
174 175
        const auto volume_f = geo_f.volume();

176 177 178
        // Get normal
        auto normal_f = entity_f.centerUnitOuterNormal();

179
        // Entity centers in references elements
180
        const auto& ref_el_i = referenceElement(geo_i);
181 182 183 184
        const auto& ref_el_o = referenceElement(geo_o);
        const auto& center_position_i  = ref_el_i.position(0,0);
        const auto& center_position_o = ref_el_o.position(0,0);
        
185 186 187
        // Get center position of the face w.r.t inside entity 
        auto face_position_i = geo_in_i.center();

188 189 190 191
        // Evaluate water content
        WaterContent water_content_i,water_content_o;
        _gf_water_content->evaluate(entity_i,center_position_i,water_content_i);
        _gf_water_content->evaluate(entity_o,center_position_o,water_content_o);
192

193 194 195 196 197 198
        // Eval flux field assume H(div) velocity field => may choose any side
        WaterFlux water_flux;
        _gf_water_flux->evaluate(entity_i, face_position_i, water_flux);
        auto water_flux_n = water_flux*normal_f;
        // auto water_vel_i = water_flux_n*water_content_i;
        // auto water_vel_o = water_flux_n*water_content_o;
199

200 201 202 203 204 205
        // Compute the effective hydrodynamic dispersion

        // bind parameterization and retrieve functions
        // _param->bind(entity_i);
        // const auto hydr_disp_func_i = _param->hydrodynamic_dispersion_f();
        // auto hydr_disp_i = hydr_disp_func_i(water_vel_i,water_content_i);
206

207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
        // _param->bind(entity_o);
        // const auto hydr_disp_func_o = _param->hydrodynamic_dispersion_f();
        // auto hydr_disp_o = hydr_disp_func_o(water_vel_o,water_content_o);

        // Evaluate diffusion coeff. from both sides
        // GradientU diff_coeff_i, diff_coeff_o;
        // hydr_disp_i.mv(normal_f,diff_coeff_i);
        // hydr_disp_o.mv(normal_f,diff_coeff_o);

        auto diff_coeff_i = _diff_coeff * water_content_i; // FIXME
        auto diff_coeff_o = _diff_coeff * water_content_o; // FIXME
        
        // Take harmonic average
        auto diff_coeff_f = 2.0/(   1.0/(diff_coeff_i + 1E-30) 
                                  + 1.0/(diff_coeff_o + 1E-30) );
222

223
        // Inside/outside solute value
224 225
        RangeU u_i = x_i(lfsu_i,0);
        RangeU u_o = x_o(lfsu_o,0);
Santiago Ospina's avatar
Santiago Ospina committed
226
        
227
        // Upwinding
228
        const auto& u_upwind = (water_flux_n>=0) ? u_i : u_o;
229 230

        // Entity centers in global coordinates
231 232
        auto center_position_i_g = geo_i.global(center_position_i);
        auto center_position_o_g = geo_o.global(center_position_o);
233 234

        // Distance between the two entity centers
235 236
        center_position_i_g -= center_position_o_g;
        auto distance = center_position_i_g.two_norm();
237 238

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

241
        // Solute flux in normal direction w.r.t the intersection
242
        auto soulte_flux_n = u_upwind*water_flux_n - diff_coeff_f*dudn;
243

244
        // Symmetric contribution to residual on inside element
245 246
        r_i.accumulate(lfsv_i, 0,  soulte_flux_n*volume_f );
        r_o.accumulate(lfsv_o, 0, -soulte_flux_n*volume_f );
247 248 249 250
    }


    /*---------------------------------------------------------------------*//**
Santiago Ospina's avatar
Santiago Ospina committed
251
     * @brief      Boundary integral depending on test and ansatz functions.
252 253 254
     *
     * @param[in]  ig      The intersection entity of the grid (inside + outside
     *                     entities)
Santiago Ospina's avatar
Santiago Ospina committed
255 256
     * @param[in]  lfsu_i  The inside ansatz local function space
     * @param[in]  x_i     The coefficients of the lfsu_i
257 258 259 260
     * @param[in]  lfsv_i  The inside test local function space
     * @param      r_i     The view of the residual vector w.r.t lfsu_i
     *
     * @tparam     IG      The type for ig
Santiago Ospina's avatar
Santiago Ospina committed
261 262 263
     * @tparam     LFSU    The type for lfsu_i
     * @tparam     X       The type for x_i
     * @tparam     LFSV    The type for lfsv_i
264 265
     * @tparam     R       The type for r_i
     */
Santiago Ospina's avatar
Santiago Ospina committed
266 267 268
    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
269
    {
270 271 272
        // Get entity dim
        // const int dim = IG::Entity::dimension;

273 274
        // Get local basis traits from local function space
        using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
275 276 277
                            Traits::LocalBasisType::Traits;

        // Get range, range field and gradient for trial space
278
        using RangeU = typename LocalBasisTraitsU::RangeType;
279 280
        // using RangeFieldU = typename LocalBasisTraitsU::RangeFieldType;
        // using GradientU = Dune::FieldVector<RangeFieldU,dim>;
281 282 283 284 285

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

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

290
        // Get geometries
291 292 293
        auto geo_f = entity_f.geometry();
        const auto& ref_el_f = referenceElement(geo_f);
        const auto& center_position_f = ref_el_f.position(0,0);
294

295 296 297 298
        // Get normal 
        auto normal_f = entity_f.centerUnitOuterNormal();

        // Get boundary condition type
299
        auto bc = _boundary->bc(entity_f,center_position_f,_time);
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
        // Get geometry of intersection in local coord. of neighbor entity
        auto geo_in_i = ig.geometryInInside();
        auto face_position_i = geo_in_i.center();

        // Evaluate convective term (i.e water flux)
        WaterFlux water_flux;
        _gf_water_flux->evaluate(entity_i, face_position_i, water_flux);
        auto water_flux_n = water_flux*normal_f;
        // auto water_vel_i = water_flux_n*water_content_i;

        // bind parameterization and retrieve functions
        // _param->bind(entity_i);
        // const auto hydr_disp_func_i = _param->hydrodynamic_dispersion_f();

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

        if (BoundaryCondition::isOutflow(bc)) 
        {
            // Check that water is going out
            if (Dune::FloatCmp::lt(water_flux_n,0.))
                return;

            // Inside unknown value
            RangeU u = x_i(lfsu_i,0);

            // Eval flux boundary condition
            auto o = _boundary->o(entity_f,center_position_f,_time);

            // convection term
            auto conv_num_flux = (u * water_flux_n + o) * volume_f;

            // integrate o
            r_i.accumulate(lfsu_i,0,conv_num_flux);
        }
        else if (BoundaryCondition::isDirichlet(bc))
337
        {
338 339 340 341 342
            // Only allow dirichlet values for influx cases
            if (Dune::FloatCmp::gt(water_flux_n,0.))
                return;

            // Get center of entity in global coordinates
343
            auto geo_i = entity_i.geometry();
344
            auto center_position_i_g = geo_i.center();
345 346

            // Face center in global coordinates
347
            auto center_position_f_g = geo_f.center();
348 349

            // Compute distance of these two points
350 351
            center_position_i_g -= center_position_f_g;
            auto distance = center_position_i_g.two_norm();
352 353

            // Get geometry of intersection in local coord. of neighbor entity
354
            auto geo_in_i = ig.geometryInInside();
355
            auto face_position_i = geo_in_i.center();
356

357 358 359
            // Evaluate water content
            WaterContent water_content_i;
            _gf_water_content->evaluate(entity_i, face_position_i, water_content_i);
360

361 362 363 364 365 366 367 368 369 370 371
            // Evaluate Dirichlet condition
            auto g = _boundary->g(entity_f,center_position_f,_time);

            // Compute the effective hydrodynamic dispersion
            // auto hydr_disp_i = hydr_disp_func_i(water_vel_i,water_content_i);

            // Calculate diffusion coefficient from inside
            // GradientU diff_coeff_i;
            // hydr_disp_i.mv(normal_f,diff_coeff_i);

            auto diff_coeff_i = _diff_coeff; // FIXME
372

Santiago Ospina's avatar
Santiago Ospina committed
373
            // Inside unknown value
374
            RangeU u = x_i(lfsu_i,0);
Santiago Ospina's avatar
Santiago Ospina committed
375

376 377 378
            // Solute flux in normal direction w.r.t the intersection
            auto soulte_flux_n = (g*water_flux_n)
                                -(diff_coeff_i*water_content_i*(g-u)/distance);
Santiago Ospina's avatar
Santiago Ospina committed
379

380 381
            // Contribution to residual from Dirichlet boundary
            r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f);
382

383 384 385
        }
        else if (BoundaryCondition::isNeumann(bc))
        {
386
            // Solute flux in normal direction w.r.t the intersection
387
            auto soulte_flux_n = _boundary->j(entity_f,center_position_f,_time);
388

389
            r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f );
390 391 392 393 394 395 396 397 398 399 400 401 402 403
        }
    }

    /*---------------------------------------------------------------------*//**
     * @brief      Sets the time.
     *
     * @param[in]  t           time of type RangeField
     *
     * @tparam     RangeField  type of the range field
     */
    template<class RangeField>
    void setTime (RangeField t)
    {
        _time = t;
404
        _gf_water_content->setTime(t);
405
        _gf_water_flux->setTime(t);
406 407
    }

408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427
    /*---------------------------------------------------------------------*//**
     * @brief      Sets the the water content grid function.
     *
     * @param[in]  gf_water_content  The water content grid function
     */
    void set_water_content(std::shared_ptr<GFWaterContent> gf_water_content)
    {
        _gf_water_content = gf_water_content;
    }

    /*---------------------------------------------------------------------*//**
     * @brief      Sets the the water flux grid function.
     *
     * @param[in]  gf_water_content  The water flux grid function
     */
    void set_water_flux(std::shared_ptr<GFWaterFlux> gf_water_flux)
    {
        _gf_water_flux = gf_water_flux;
    }

428
private:
429
    const std::shared_ptr<const Boundary> _boundary;
430 431
    std::shared_ptr<GFWaterFlux> _gf_water_flux;
    std::shared_ptr<GFWaterContent> _gf_water_content;
432
    double _time;
433
    double _diff_coeff;
434 435 436 437 438
};

/**
 * @brief      Temporal local operator for the transport equation in unsaturated
 *             media in a finite volume scheme.
439
 * @details    It solves the temporal part of the transport equation:
Lukas Riedel's avatar
Lukas Riedel committed
440 441
 *
 *             @f{eqnarray*}{
442 443
 *             \partial_t[\theta C_w] +
 *             \nabla\cdot [\textbf{j}_w C_w] -
Lukas Riedel's avatar
Lukas Riedel committed
444
 *             \nabla [\theta \mathsf{D}_{\text{eff}}\nabla C_w]&=0 &\qquad \text{in }
445
 *             \Omega\\
Lukas Riedel's avatar
Lukas Riedel committed
446
 *             C_w &= g &\qquad \text{on } \Gamma_D
447
 *             \subseteq\partial\Omega\\
Lukas Riedel's avatar
Lukas Riedel committed
448
 *             \nabla C_w \cdot \textbf{n} &= \textbf{j}_{\scriptscriptstyle
449
 *             C_w}& \qquad \text{on } \Gamma_N =\partial\Omega \backslash
Lukas Riedel's avatar
Lukas Riedel committed
450 451 452
 *             \Gamma_D
 *             @f}
 *
453 454
 * @author     Santiago Ospina De Los Ríos
 * @date       2018
455
 * @ingroup    LocalOperators
456
 * @ingroup    TransportModel
457 458
 *
 * @bug        The water content is not being updated by the set time to the
459
 *             right state of the grid function.
460
 *
461 462
 * @tparam     GFWaterContent  Type of a grid function which provides the water
 *                             content
463
 */
464
template<class GFWaterContent>
465 466
  class TransportFVTemporalOperator
  : public Dune::PDELab::NumericalJacobianVolume<
467
                            TransportFVTemporalOperator<GFWaterContent>>
468
  , public Dune::PDELab::NumericalJacobianApplyVolume<
469
                            TransportFVTemporalOperator<GFWaterContent>>
470 471 472 473 474 475 476 477 478 479 480 481 482 483 484
  , public Dune::PDELab::FullSkeletonPattern
  , public Dune::PDELab::FullVolumePattern
  , public Dune::PDELab::LocalOperatorDefaultFlags
  , public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double>
{
public:

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

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

public:

485 486 487
    TransportFVTemporalOperator()
    : _time(0.)
    { }
488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507

    /*---------------------------------------------------------------------*//**
     * @brief      Volume integral depending on test and ansatz functions
     *
     * @param[in]  eg    THe entity of the grid
     * @param[in]  lfsu  The ansatz local function space
     * @param[in]  x     The coefficients of the lfsu
     * @param[in]  lfsv  The test local function space
     * @param      r     The view of the residual vector w.r.t lfsu
     *
     * @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
     */
    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
    {
508 509 510 511 512
        // 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;
Santiago Ospina's avatar
Santiago Ospina committed
513

514
        RangeU u = x(lfsu,0);
Santiago Ospina's avatar
Santiago Ospina committed
515 516 517 518

        // entity geometry
        auto geo = eg.geometry();

519
        // inside cell center
520
        const auto& ref_el = referenceElement(geo);
521
        const auto& center_position = ref_el.position(0,0);
522

523
        typename GFWaterContent::Traits::RangeType water_content;
524

525
        _gf_water_content->evaluate(eg.entity(),center_position,water_content);
526 527

        // update residual
528
        r.accumulate(lfsv ,0 , water_content*u*geo.volume());
529 530 531 532 533 534 535 536 537 538 539 540 541
    }

    /*---------------------------------------------------------------------*//**
     * @brief      Sets the time.
     *
     * @param[in]  t           time of type RangeField
     *
     * @tparam     RangeField  type of the range field
     */
    template<class RF>
    void setTime (RF t)
    {
        _time = t;
542
        _gf_water_content->setTime(t);
543 544
    }

545 546 547 548 549 550 551 552 553 554
    /*---------------------------------------------------------------------*//**
     * @brief      Sets the the water content grid function.
     *
     * @param[in]  gf_water_content  The water content grid function
     */
    void set_water_content(std::shared_ptr<GFWaterContent> gf_water_content)
    {
        _gf_water_content = gf_water_content;
    }

555
private:
556
    std::shared_ptr<GFWaterContent> _gf_water_content;
557 558 559 560 561 562 563 564 565
    double _time;
};

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


#endif // DUNE_DORIE_TRANSPORT_OPERATOR_HH