test-mass-conservation.hh 9.96 KB
Newer Older
1 2 3 4
#ifndef TEST_MASS_CONSERVATION_HH
#define TEST_MASS_CONSERVATION_HH

#include <algorithm>
5 6
#include <vector>
#include <memory>
7 8
#include <iostream>
#include <algorithm>
9

10 11
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/exceptions.hh>
12 13
#include <dune/common/parametertree.hh>
#include <dune/common/parallel/mpihelper.hh>
14 15

#include <dune/geometry/referenceelements.hh>
16

17
#include <dune/pdelab/common/quadraturerules.hh>
18 19
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/gridoperator/gridoperator.hh>
20
#include <dune/pdelab/constraints/common/constraintstransformation.hh>
21 22 23 24 25 26
#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>

27
#include <dune/dorie/common/typedefs.hh>
28
#include <dune/dorie/common/logging.hh>
29
#include <dune/dorie/common/boundary_condition/neumann.hh>
30
#include <dune/dorie/model/richards/richards.cc>
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76

namespace Dune {
namespace Dorie {

/// Local Operator evaluating the water content of a grid cell
template <typename Traits, typename Parameter>
class WaterContentOperator : public Dune::PDELab::FullVolumePattern,
                             public Dune::PDELab::LocalOperatorDefaultFlags
{
public:
  // residual assembly flags
  enum { doAlphaVolume = true };

private:
  static constexpr int dim = Traits::dim;
  using RF = typename Traits::RF;
  using Scalar = typename Traits::Scalar;

  const Parameter &param;      //!< class containing parameterization
  const int quadrature_factor; //!< factor to FEM integration order
  const int intorderadd;       //!< integration order addend

public:
  WaterContentOperator(
    const Parameter &_param,
    const int _quadrature_factor = 2,
    const int _intorderadd = 2)
    : param(_param),
      quadrature_factor(_quadrature_factor),
      intorderadd(_intorderadd)
  { }

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

77 78 79 80 81
    // 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();

82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
    // loop over quadrature points
    for (const auto &it : quadratureRule(gt, intorder))
    {
      // evaluate position in element local and global coordinates
      const auto p_local = it.position();

      // evaluate basis function
      std::vector<Scalar> phi(lfsu.size());
      lfsu.finiteElement().localBasis().evaluateFunction(p_local, phi);

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

      // calculate water content from matric head
98
      const RF water_content = water_content_f(saturation_f(u));
99 100 101 102 103 104 105 106 107 108 109 110 111 112 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 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158

      // integration factor
      const RF factor = it.weight() * gt.integrationElement(p_local);

      // update residual
      r.accumulate(lfsv, 0, water_content * factor);
    } // quadrature rule
  }   // alpha_volume
};

/// Local Operator evaluating the boundary flux for each cell
/** This only includes fluxed through Neumann boundary conditions
 */
template <typename Traits, typename Parameter, typename Boundary>
class FluxOperator : public Dune::PDELab::FullVolumePattern,
                     public Dune::PDELab::LocalOperatorDefaultFlags,
                     public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<typename Traits::RangeField>
{
public:
  // residual assembly flags
  enum { doLambdaBoundary = true };

private:
  static constexpr int dim = Traits::dim;
  using RF = typename Traits::RF;
  using Scalar = typename Traits::Scalar;

  const Parameter &param;      //!< class containing parameterization
  const Boundary &boundary;    //!< class containing BC information
  const int quadrature_factor; //!< factor to FEM integration order
  const int intorderadd;       //!< integration order addend

  RF time;

public:
  FluxOperator(
    const Parameter &_param,
    const Boundary &_boundary,
    const int _quadrature_factor = 2,
    const int _intorderadd = 2)
  : param(_param),
    boundary(_boundary),
    quadrature_factor(_quadrature_factor),
    intorderadd(_intorderadd)
  { }

  void setTime(RF t) { time = t; }

  template <typename IG, typename LFSV, typename R>
  void lambda_boundary(
      const IG &ig,
      const LFSV &lfsv_s,
      R &r_s) const
  {
    const int order = lfsv_s.finiteElement().localBasis().order();
    const int intorder = intorderadd + quadrature_factor * order;

    // get element geomentry
    auto gtface = ig.geometry();

159 160 161 162 163 164 165 166 167 168 169 170 171 172
    // check boundary condition type
    const auto bc = boundary.bc(ig.intersection());
    if (bc->type() != Operator::BCType::Neumann)
      return;

    // evaluate flux boundary condition
    RF normal_flux = bc->evaluate(time);

    // adjust flux to surface tilt
    const auto& bc_neumann
      = dynamic_cast<const NeumannBoundaryCondition<RF>&>(*bc);
    if (bc_neumann.horizontal_projection())
      normal_flux *= std::abs( ig.centerUnitOuterNormal() * param.gravity() );

173 174 175 176 177 178 179
    // loop over quadrature points
    for (const auto &it : quadratureRule(gtface, intorder))
    {
      // integration factor
      const RF factor = it.weight() * ig.geometry().integrationElement(it.position());

      // update residual
180
      r_s.accumulate(lfsv_s, 0, - normal_flux * factor);
181 182 183 184
    } // quadrature rule
  }   // lambda_boundary
};

185
/// Derived Model for testing mass conservation. Assumes
186
template <typename Traits>
187
class ModelTest : public ModelRichards<Traits>
188 189 190
{
protected:
  // fetch typedefs from base class
191
  using Sim = ModelRichards<Traits>;
192 193 194 195 196 197
  using GV = typename Traits::GV;
  using RF = typename Traits::RF;
  using GFS = typename Traits::GFS;
  using MBE = typename Traits::MBE;
  using Grid = typename Traits::Grid;
  using U = typename Traits::U;
198

199
  using GridCreator = typename Traits::GridCreator;
200

201
  using FlowParameters = typename Traits::FlowParameters;
202
  using FlowBoundary = typename Traits::FlowBoundary;
203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219

  /// Error Function Space Helper
  using ERRGFSHelper = Dune::Dorie::LinearSolverGridFunctionSpaceHelper<GV, RF, Traits::GridGeometryType>;
  /// Linear solver GFS
  using ERRGFS = typename ERRGFSHelper::Type;
  /// LSGFS Constraints Container
  using ERRCC = typename ERRGFSHelper::CC;
  /// Solution vector type
  using U0 = Dune::PDELab::Backend::Vector<ERRGFS, RF>;

  /// Empty constraints container
  using NoTrafo = Dune::PDELab::EmptyTransformation;
  /// Grid operator for Error LOP
  template <typename LOP>
  using ERRGO = Dune::PDELab::GridOperator<GFS, ERRGFS, LOP, MBE, RF, RF, RF, NoTrafo, NoTrafo>;

public:
220
  /// Construct this model instance.
221 222
  /** \todo Give error limits here
   */
223
  ModelTest(
224
    Dune::ParameterTree& _inifile,
225
		const GridCreator& _grid_creator,
226
		Dune::MPIHelper& _helper
227
  ) : ModelRichards<Traits>(_inifile, _grid_creator, _helper)
228 229
    , acc(0.)
    , acc_square(0.)
230 231
  { 
    // disable output
232
    // this->set_policy(OutputPolicy::None);
233
  }
234 235 236 237

  /// Compute the water flux into a certain domain
  double water_flux(U &solution, const RF time)
  {
238 239 240
    using FluxLOP = typename Dune::Dorie::FluxOperator<Traits,
                                                       FlowParameters, 
                                                       FlowBoundary>;
241 242 243 244
    using FluxGO = ERRGO<FluxLOP>;

    // create GFS and operators
    auto errgfs = ERRGFSHelper::create(this->gv);
245
    FluxLOP fluxlop(*this->fparam, *this->fboundary);
246 247 248 249 250 251 252 253 254 255 256 257 258 259
    fluxlop.setTime(time);
    MBE mbe(0);
    FluxGO fluxgo(*this->gfs, errgfs, fluxlop, mbe);
    U0 eta(errgfs, 0.0);

    // assemble residual and accumulate value
    fluxgo.residual(solution, eta);
    auto &eta_n = Dune::PDELab::Backend::native(eta);
    return std::accumulate(eta_n.begin(), eta_n.end(), 0.0);
  }

  /// Compute the total water content of a certain solution
  double water_content(U &solution)
  {
260 261
    using WCLOP = typename Dune::Dorie::WaterContentOperator<Traits, 
                                                             FlowParameters>;
262 263 264 265
    using WCGO = ERRGO<WCLOP>;

    // create GFS and operators
    auto errgfs = ERRGFSHelper::create(this->gv);
266
    WCLOP wclop(*this->fparam);
267 268 269 270 271 272 273 274 275 276
    MBE mbe(0);
    WCGO wcgo(*this->gfs, errgfs, wclop, mbe);
    U0 eta(errgfs, 0.0);

    // assemble residual and accumulate value
    wcgo.residual(solution, eta);
    auto &eta_n = Dune::PDELab::Backend::native(eta);
    return std::accumulate(eta_n.begin(), eta_n.end(), 0.0);
  }

277
  void step() override
278
  {
279 280
    const auto time_old = this->_time;
    const double wc_old = water_content(*this->u);
281
    // compute expected integrated flux (water content change)
282
    const double flux = water_flux(*this->u, time_old);
283

284
    ModelRichards<Traits>::step();
285

286 287 288 289
    const auto time_new = this->_time;
    const auto dt = time_new - time_old;
    const auto flux_int = flux * dt;

290
    // compute new water content
291
    const auto wc_new = water_content(*this->u);
292 293 294 295 296 297 298 299 300
    const auto deviation = (wc_new - wc_old) - flux_int;
    acc += deviation;
    acc_square += deviation * deviation;

    std::cout << "wc_old: " << wc_old << std::endl;
    std::cout << "wc_new: " << wc_new << std::endl;
    std::cout << "integrated flux: " << flux_int << std::endl;
    std::cout << "deviation: " << deviation << std::endl;
  }
301

302 303 304 305
  double run_test ()
  {
    acc = 0.0;        // accumulated deviation
    acc_square = 0.0; // accumulated squared deviation
306

307
    this->run();
308 309 310 311 312 313

    std::cout << "total deviation: " << acc << std::endl;
    std::cout << "total normalized deviation: " << std::sqrt(acc_square) << std::endl;

    return acc;
  }
314 315 316 317

private:
  double acc;        // accumulated deviation
  double acc_square; // accumulated squared deviation
318 319 320 321 322
};

} // namespace Dorie
} // namespace Dune

323
#endif // TEST_MASS_CONSERVATION_HH