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:
25 26 27 28 29 30 31 32 33 34 35
 * @f{eqnarray*}{
 *             \partial_t[\theta C_w] +
 *             \nabla\cdot [\textbf{j}_w C_w] -
 *             \nabla [\theta \mathsf{D}_{eff}\nabla C_w]=0 &\qquad \text{in }
 *             \Omega\\ 
 *             C_w = g &\qquad \text{on } \Gamma_D
 *             \subseteq\partial\Omega\\
 *             \nabla C_w \cdot \textbf{n} = \textbf{j}_{\scriptscriptstyle
 *             C_w}& \qquad \text{on } \Gamma_N =\partial\Omega \backslash
 *             \Gamma_D
 * @f}
36 37 38 39
 * @author     Santiago Ospina De Los Ríos
 * @date       2018
 * @copyright  MIT License.
 * @ingroup    LocalOperators
40
 * 
41 42
 * @todo        Use diffusion coefficient grid function.
 * @todo        Implement outflow boundary condition.
43
 * 
44
 *
45 46 47 48 49 50
 * @tparam     Boundary                  Type of the class providing boundary
 *                                       conditions
 * @tparam     GridFunctionWaterFlux     Type of a grid function which provides
 *                                       the water flux
 * @tparam     GridFunctionWaterContent  Type of a grid function which provides
 *                                       the water content
51
 */
52
template<class Boundary, class GridFunctionWaterFlux, class GridFunctionWaterContent>
53 54 55 56 57
class TransportFVSpatialOperator
    : public Dune::PDELab::NumericalJacobianVolume<
                        TransportFVSpatialOperator<
                            Boundary,
                            GridFunctionWaterFlux,
58
                            GridFunctionWaterContent
59
                        >
60
                    > 
61 62 63 64
    , public Dune::PDELab::NumericalJacobianSkeleton<
                        TransportFVSpatialOperator<
                            Boundary,
                            GridFunctionWaterFlux,
65
                            GridFunctionWaterContent
66 67 68 69 70 71
                        >
                    >
    , public Dune::PDELab::NumericalJacobianBoundary<
                        TransportFVSpatialOperator<
                            Boundary,
                            GridFunctionWaterFlux,
72
                            GridFunctionWaterContent
73 74 75 76 77 78 79 80 81 82 83 84 85
                        >
                    >
    , 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
86
    // enum { doAlphaVolume    = true };
87 88
    enum { doAlphaBoundary = true };
    enum { doAlphaSkeleton  = true };
89 90
    // enum { doLambdaVolume   = true };
    // enum { doLambdaBoundary   = true };
91
private:
92 93
    static_assert(std::is_same<
            typename GridFunctionWaterFlux::Traits::GridViewType,
94
            typename GridFunctionWaterContent::Traits::GridViewType>::value, 
95
        "TransportFVSpatialOperator: GridFunctionWaterFlux and"
96
        "GridFunctionWaterContent has to use the same grid view.");
97

98
    using WaterFlux =typename GridFunctionWaterFlux::Traits::RangeType;
99
    using WaterContent =typename GridFunctionWaterContent::Traits::RangeType;
100 101

public:
102

103 104 105
    /**
     * @brief      Constructor of TransportFVSpatialOperator
     *
106 107 108 109
     * @param      boundary          The boundary terms
     * @param[in]  gf_water_flux     The grid function of water flux
     * @param[in]  gf_water_content  The grid function of water content
     * @param[in]  diff_coeff        The diffusion coefficient
110
     */
111 112
    TransportFVSpatialOperator(
        Boundary& boundary, 
113 114
        std::shared_ptr<GridFunctionWaterFlux> gf_water_flux, 
        std::shared_ptr<GridFunctionWaterContent> gf_water_content,
115
        double diff_coeff
116
    )   : _boundary(boundary)
117 118
        , _gf_water_flux(gf_water_flux)
        , _gf_water_content(gf_water_content)
119
        , _diff_coeff(diff_coeff)
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
    {}

    /*---------------------------------------------------------------------*//**
     * @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>
138 139
    void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, 
                                     const LFSV& lfsv, R& r) const
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 167 168
    {}

    /*---------------------------------------------------------------------*//**
     * @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
    {
169 170 171 172 173 174 175 176 177 178 179 180 181 182
        // 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;

        // 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();
183 184
        const auto& entity_i = ig.inside();
        const auto& entity_o = ig.outside();
185 186

        // Get geometries
187
        auto geo_f = entity_f.geometry();
188 189
        auto geo_i = entity_i.geometry();
        auto geo_o = entity_o.geometry();
190

191
        // Get geometry of intersection in local coord. of neighbor entities
192
        auto geo_in_i = entity_f.geometryInInside();
193

194 195 196 197 198 199 200 201 202 203 204 205 206
        // Get volume of entities
        const auto volume_f = geo_f.volume();

        // Entity centers in references elements
        const auto& ref_el_i  = referenceElement(geo_i);
        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);
        
        // 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);
207

208
        // Evaluate diffusion coeff. from both sides and take harmonic average
209 210 211 212 213
        auto diff_coeff_i = _diff_coeff * water_content_i; // FIXME: evaluate from grid function
        auto diff_coeff_o = _diff_coeff * water_content_o; // FIXME: evaluate from grid function
        auto normal_f = entity_f.centerUnitOuterNormal();
        auto diff_coeff_f = 2.0/(   1.0/(diff_coeff_i + 1E-30) 
                                  + 1.0/(diff_coeff_o + 1E-30) );
214 215

        // Get center position of the face w.r.t inside entity 
216
        auto face_position_i = geo_in_i.center();
217

218
        // Evaluate velocity field assume H(div) velocity field => may choose any side
219
        WaterFlux water_flux_i;
220 221
        _gf_water_flux->evaluate(entity_i, face_position_i, water_flux_i);
        auto water_flux_n = water_flux_i*normal_f;
222

223
        // Inside/outside solute value
224 225
        RangeU u_i = x_i(lfsu_i,0);
        RangeU u_o = x_o(lfsu_o,0);
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
    }


    /*---------------------------------------------------------------------*//**
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)
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
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
     */
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 273 274 275 276 277 278 279
        // 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;

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

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

284
        // Get geometries
285 286 287
        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);
288

289
        auto bc = _boundary.bc(entity_f,center_position_f,_time);
290 291 292

        if (BoundaryCondition::isDirichlet(bc))
        {
293
            auto geo_i = entity_i.geometry();
294
            auto center_position_i_g = geo_i.center();
295 296

            // Face center in global coordinates
297
            auto center_position_f_g = geo_f.center();
298 299

            // Compute distance of these two points
300 301
            center_position_i_g -= center_position_f_g;
            auto distance = center_position_i_g.two_norm();
302 303

            // Evaluate Dirichlet condition
304
            auto g = _boundary.g(entity_f,center_position_f,_time);
305 306

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

            // Get geometry of intersection in local coord. of neighbor entity
310
            auto geo_in_i = ig.geometryInInside();
311
            auto face_position_i = geo_in_i.center();
312 313

            // Evaluate convective term (i.e water flux)
314
            WaterFlux water_flux_i;
315 316 317
            _gf_water_flux->evaluate(entity_i, face_position_i, water_flux_i);
            auto normal_f = ig.unitOuterNormal(center_position_f);
            auto water_flux_n = water_flux_i*normal_f;
318

319 320 321
            // Evaluate water content
            WaterContent water_content_i;
            _gf_water_content->evaluate(entity_i, face_position_i, water_content_i);
322 323

            // Evaluate diffusion coefficient from inside
324
            auto diff_coeff_i = _diff_coeff; // FIXME: evaluate from grid function
325

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

Santiago Ospina's avatar
Santiago Ospina committed
329 330 331
            // Only allow dirichlet values for influx cases
            if (water_flux_n<=0){
                // Solute flux in normal direction w.r.t the intersection
332
                auto soulte_flux_n = (g*water_flux_n)-(diff_coeff_i*water_content_i*(g-u)/distance);
Santiago Ospina's avatar
Santiago Ospina committed
333 334

                // Contribution to residual from Dirichlet boundary
335
                r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f);
Santiago Ospina's avatar
Santiago Ospina committed
336 337
            }

338

339 340 341 342
        }
        else if (BoundaryCondition::isNeumann(bc))
        {
            // Get geometry of intersection in local coord. of neighbor entity
343
            auto geo_in_i = ig.geometryInInside();
344
            auto face_position_i = geo_in_i.center();
345

346 347 348
            // Evaluate water content
            typename GridFunctionWaterContent::Traits::RangeType water_content_i;
            _gf_water_content->evaluate(entity_i, face_position_i, water_content_i);
349 350

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

            // Evaluate diffusion coefficient from inside
354
            auto diff_coeff_i = _diff_coeff;
355 356

            // Contribution to residual from Neumann boundary
357
            auto j = _boundary.j(entity_f,center_position_f,_time);
358 359

            // Solute flux in normal direction w.r.t the intersection
360
            auto soulte_flux_n = -diff_coeff_i*j*water_content_i;
361

362
            r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f );
363 364 365
        }
    }

366 367 368 369 370 371 372 373 374 375 376 377 378
    /*---------------------------------------------------------------------*//**
     * @brief      Residual contribution of volume integral
     *
     * @param[in]  eg    THe entity of the grid
     * @param[in]  lfsu  The ansatz local function space
     * @param      r     The view of the residual vector w.r.t lfsu
     * 
     * @tparam     EG      The type for eg
     * @tparam     LFSV    The type for lfsv and
     * @tparam     R       The type for r
     */
    template<typename EG, typename LFSV, typename R>
    void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
Santiago Ospina's avatar
Santiago Ospina committed
379
    {}
380 381

    /*---------------------------------------------------------------------*//**
382
     * @brief      Residual contribution of boundary integral
383 384 385 386 387 388 389
     *
     * @param[in]  ig      The intersection entity of the grid (inside + outside
     *                     entities)
     * @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
390
     * @tparam     LFSV    The type for lfsv_i and
391 392
     * @tparam     R       The type for r_i
     */
393 394 395 396 397 398 399 400 401 402 403 404 405
    template<typename IG, typename LFSV, typename R>
    void lambda_boundary (const IG& ig, const LFSV& lfsv_i, R& r_i) const
    {}

    /*---------------------------------------------------------------------*//**
     * @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)
406
    {
407 408 409
        _time = t;
    }

410 411
private:
    Boundary& _boundary;
412 413
    const std::shared_ptr<GridFunctionWaterFlux> _gf_water_flux;
    const std::shared_ptr<GridFunctionWaterContent> _gf_water_content;
414
    double _time;
415
    double _diff_coeff;
416 417 418 419 420
};

/**
 * @brief      Temporal local operator for the transport equation in unsaturated
 *             media in a finite volume scheme.
421 422 423 424 425 426 427 428 429 430 431 432 433 434 435
 * @details    It solves the temporal part of the transport equation:
 * @f{eqnarray*}{
 *             \partial_t[\theta C_w] +
 *             \nabla\cdot [\textbf{j}_w C_w] -
 *             \nabla [\theta \mathsf{D}_{eff}\nabla C_w]=0 &\qquad \text{in }
 *             \Omega\\ 
 *             C_w = g &\qquad \text{on } \Gamma_D
 *             \subseteq\partial\Omega\\
 *             \nabla C_w \cdot \textbf{n} = \textbf{j}_{\scriptscriptstyle
 *             C_w}& \qquad \text{on } \Gamma_N =\partial\Omega \backslash
 *             \Gamma_D
 * @f}
 * @author     Santiago Ospina De Los Ríos
 * @date       2018
 * @copyright  MIT License.
436
 * @ingroup    LocalOperators
437 438 439
 * 
 * @bug        The water content is not being updated by the set time to the 
 *             right state of the grid function.
440 441 442
 *
 * @tparam     GridFunctionWaterContent  Type of a grid function which provides
 *                                       the water content
443
 */
444
template<class GridFunctionWaterContent>
445 446
  class TransportFVTemporalOperator
  : public Dune::PDELab::NumericalJacobianVolume<
447
                            TransportFVTemporalOperator<GridFunctionWaterContent>>
448
  , public Dune::PDELab::NumericalJacobianApplyVolume<
449
                            TransportFVTemporalOperator<GridFunctionWaterContent>>
450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465
  , 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:

    TransportFVTemporalOperator(
466 467
        std::shared_ptr<GridFunctionWaterContent> gf_water_content
    )   : _gf_water_content(gf_water_content)
468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488
    {}

    /*---------------------------------------------------------------------*//**
     * @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
    {
489 490 491 492 493
        // 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;
494

495
        RangeU u = x(lfsu,0);
496 497 498 499

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

500
        // inside cell center
501
        const auto& ref_el = referenceElement(geo);
502
        const auto& center_position = ref_el.position(0,0);
503

504
        typename GridFunctionWaterContent::Traits::RangeType water_content;
505

506
        _gf_water_content->evaluate(eg.entity(),center_position,water_content);
507 508

        // update residual
509
        r.accumulate(lfsv ,0 , water_content*u*geo.volume());
510 511 512 513 514 515 516 517 518 519 520 521 522
    }

    /*---------------------------------------------------------------------*//**
     * @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;
523
        _gf_water_content->setTime(t);
524 525 526
    }

private:
527
    std::shared_ptr<GridFunctionWaterContent> _gf_water_content;
528 529 530 531 532 533 534 535 536
    double _time;
};

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


#endif // DUNE_DORIE_TRANSPORT_OPERATOR_HH