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.)

transport.hh 24 KB
Newer Older
1 2
#ifndef DUNE_DORIE_MODEL_TRANSPORT_HH
#define DUNE_DORIE_MODEL_TRANSPORT_HH
3 4 5 6 7

#include <memory>
#include <string>
#include <iostream>
#include <limits>
8 9
#include <algorithm>
#include <utility>
10 11 12

#include <dune/common/parametertree.hh>
#include <dune/common/parallel/mpihelper.hh>
13
#include <dune/common/float_cmp.hh>
14 15 16 17 18 19

#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/gridoperator/gridoperator.hh>
#include <dune/pdelab/gridoperator/onestep.hh>
#include <dune/pdelab/stationary/linearproblem.hh>

20 21 22
#include <dune/dorie/common/typedefs.hh>
#include <dune/dorie/common/utility.hh>
#include <dune/dorie/common/gfs_helper.hh>
23
#include <dune/dorie/common/interpolator.hh>
24
#include <dune/dorie/common/time_step_controller.hh>
25
#include <dune/dorie/common/grid_creator.hh>
26 27
#include <dune/dorie/common/grid_function_writer.hh>
#include <dune/dorie/common/cfl_condition.hh>
28
#include <dune/dorie/common/flux_reconstruction/raviart_thomas.hh>
29
#include <dune/dorie/common/boundary_condition/manager.hh>
30
#include <dune/dorie/common/logging.hh>
31

32
#include <dune/dorie/model/base.hh>
33
#include <dune/dorie/model/transport/adapters/total_solute.hh>
34 35
#include <dune/dorie/model/transport/adapters/peclet.hh>
#include <dune/dorie/model/transport/adapters/effective_hydrodynamic_dispersion.hh>
36
#include <dune/dorie/model/transport/solute_parameters.hh>
37
#include <dune/dorie/model/transport/initial_condition.hh>
38
#include <dune/dorie/model/transport/local_operator_FV.hh>
39
#include <dune/dorie/model/transport/local_operator_DG.hh>
40

41 42 43 44
#ifdef GTEST
  #include <gtest/gtest.h>
#endif

45 46 47 48
namespace Dune{
namespace Dorie{

/*-------------------------------------------------------------------------*//**
49
 * @brief      Traits for the class ModelTransport.
50
 * @details    This class extends BaseTraits to be used for the
51
 *             implementation of the Transport model.
52 53
 * @author     Santiago Ospina De Los Ríos
 * @date       2018
54
 * @ingroup    TransportModel
55
 *
56
 * @todo       Implement higher order methods.
57
 *
58
 * @tparam     BaseTraits          Traits defining domain and range field
59
 *                                 properties of the model.
60 61
 * @tparam     GFWaterFluxType     Grid Function type of the water flux
 * @tparam     GFWaterContentType  Grid Function type of the water content
62 63 64
 * @tparam     k                   Order of the polynomial degree used for the
 *                                 basis functions in the DG method.
 *                                 Zero (0) enables the FV method.
65
 */
66 67 68 69
template<class BaseTraits,
         class GFWaterFluxType,
         class GFWaterContentType,
         unsigned int k>
70
struct ModelTransportTraits : public BaseTraits
71 72
{
  static constexpr int dim = BaseTraits::dim;
73
  static constexpr int order = k;
74
  static constexpr int flux_order = (order==0) ? 0 : (order-1);
75 76 77 78 79
  using RF = typename BaseTraits::RF;
  using Grid = typename BaseTraits::Grid;
  using GV = typename BaseTraits::GV;

  /// (Inestationary) Grid Function Type for water flux
80
  static_assert(GFWaterFluxType::Traits::dimRange == dim,
81 82 83
    "Water Flux Grid Function has wrong range");
  using GFWaterFlux = GFWaterFluxType;
  /// (Inestationary) Grid Function Type for water content
84
  static_assert(GFWaterContentType::Traits::dimRange == 1,
85 86
    "Water Content Grid Function has wrong range");
  using GFWaterContent = GFWaterContentType;
87 88

  /// GFS Helper
89 90 91 92 93 94 95 96
  using GFSHelper = std::conditional_t<
    order==0,
    Dune::Dorie::LinearSolverGridFunctionSpaceHelper<
      GV, RF, BaseTraits::GridGeometryType>,
    Dune::Dorie::GridFunctionSpaceHelper<
      GV, RF, order, BaseTraits::GridGeometryType>
  >;

97 98 99 100
  /// Problem GFS
  using GFS = typename GFSHelper::Type;
  /// GFS Constraints Container
  using CC = typename GFSHelper::CC;
101 102 103 104 105 106 107 108 109 110 111
  /// Finite Element Map Type
  using FEM = typename GFSHelper::FEM;

  /// LSGFS Helper
  using LSGFSHelper = Dune::Dorie::LinearSolverGridFunctionSpaceHelper<
    GV, RF, BaseTraits::GridGeometryType
  >;
  /// Linear solver GFS
  using LSGFS = typename LSGFSHelper::Type;
  /// LSGFS Constraints Container
  using LSCC = typename LSGFSHelper::CC;
112

113 114
  // -- DORiE Class Definitions -- //
  /// Wrapper for grid and volume & boundary mappings
115
  using GridCreator = Dune::Dorie::GridCreator<typename BaseTraits::Grid>;
116 117
  /// Class defining the solute parameterization
  using SoluteParameters = Dune::Dorie::SoluteParameters<BaseTraits>;
118
  /// Class defining boundary conditions
119 120
  using SoluteBoundary = Dune::Dorie::BoundaryConditionManager<
    BaseTraits, BCMMode::transport>;
121
  /// Class defining the initial condition
122 123
  using SoluteInitial = Dune::Dorie::InitialCondition<BaseTraits>;
  /// Class defining the initial condition factory
124 125
  using SoluteInitialFactory
    = Dune::Dorie::TransportInitialConditionFactory<BaseTraits>;
126

127
  /// Local spatial operator
128 129 130 131 132 133 134 135
  using SLOP = std::conditional_t<
    order==0,
    Dune::Dorie::Operator::TransportFVSpatialOperator<
      SoluteParameters, SoluteBoundary, GFWaterFlux, GFWaterContent>,
    Dune::Dorie::Operator::TransportDGSpatialOperator<
      FEM, SoluteParameters, SoluteBoundary, GFWaterFlux, GFWaterContent>
  >;

136
  /// Local temporal operator
137 138 139 140 141 142
  using TLOP = std::conditional_t<
    order==0,
    Dune::Dorie::Operator::TransportFVTemporalOperator<GFWaterContent>,
    Dune::Dorie::Operator::TransportDGTemporalOperator<FEM, GFWaterContent>
  >;

143
  /// Time controller
144
  using TimeStepController = Dune::Dorie::TimeStepController<double>;
145 146 147 148 149 150 151 152 153 154 155

  // -- DUNE Class Defintions -- //
  /// Matrix backend type
  using MBE = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
  /// Grid operator for spatial LOP
  using GO0 = Dune::PDELab::GridOperator<GFS,GFS,SLOP,MBE,RF,RF,RF,CC,CC>;
  /// Grid operator for temporal LOP
  using GO1 = Dune::PDELab::GridOperator<GFS,GFS,TLOP,MBE,RF,RF,RF,CC,CC>;

  // -- setup -- //
  /// One step grid operator for merging both GOs
156 157
  using ImplicitIGO = Dune::PDELab::OneStepGridOperator<GO0, GO1, true>;
  using ExplicitIGO = Dune::PDELab::OneStepGridOperator<GO0, GO1, false>;
158
  /// Solution vector type
159
  using U = typename ImplicitIGO::Traits::Domain;
160
  /// Linear solver type
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
  using ImplicitLS = std::conditional_t<
    order == 0,
    Dune::PDELab::ISTLBackend_BCGS_AMG_SSOR<ImplicitIGO>,
    Dune::PDELab::ISTLBackend_OVLP_AMG_4_DG<ImplicitIGO,
                                            CC,
                                            LSGFS,
                                            LSCC,
                                            Dune::PDELab::CG2DGProlongation,
                                            Dune::SeqSSOR,
                                            Dune::BiCGSTABSolver>>;
  using ExplicitLS = std::conditional_t<
    order == 0,
    Dune::PDELab::ISTLBackend_BCGS_AMG_SSOR<ExplicitIGO>,
    Dune::PDELab::ISTLBackend_OVLP_AMG_4_DG<ExplicitIGO,
                                            CC,
                                            LSGFS,
                                            LSCC,
                                            Dune::PDELab::CG2DGProlongation,
                                            Dune::SeqSSOR,
                                            Dune::BiCGSTABSolver>>;
181
  /// Methods computing the time step
182 183 184 185 186 187
  using ImplicitSLPS
    = Dune::PDELab::StationaryLinearProblemSolver<ImplicitIGO, ImplicitLS, U>;
  using ImplicitOSM
    = Dune::PDELab::OneStepMethod<RF, ImplicitIGO, ImplicitSLPS, U, U>;
  using ExplicitOSM
    = Dune::PDELab::ExplicitOneStepMethod<RF, ExplicitIGO, ExplicitLS, U, U>;
188

189
  /// Model state provider
190 191 192 193 194 195 196 197 198 199 200
  struct ConstState 
  {
    using GridFunctionSpace = const GFS;
    using Coefficients = const U;
    using TimeField = const typename BaseTraits::TimeField;

    std::shared_ptr<GridFunctionSpace> grid_function_space;
    std::shared_ptr<Coefficients> coefficients;
    TimeField time;
  };

201
  /// Model state provider
202 203 204 205
  struct State 
  {
    using GridFunctionSpace = GFS;
    using Coefficients = U;
206
    using TimeField = typename BaseTraits::TimeField;
207 208 209

    std::shared_ptr<GridFunctionSpace> grid_function_space;
    std::shared_ptr<Coefficients> coefficients;
210 211 212 213 214 215
    TimeField time;

    // Implicit conversion to constant state
    operator ConstState() {
      return ConstState{grid_function_space,coefficients,time};
    }
216 217
  };

218 219 220
  // -- Grid Functions of state valriables -- //
  /// Solute
  using GFSolute = Dune::PDELab::DiscreteGridFunction<GFS,U>;
221
  /// TotalSolute
222
  using GFTotalSolute = Dune::Dorie::TotalSoluteAdapter<GFSolute,GFWaterContent>;
223
  /// Water Flux reconstruction
224
  static constexpr bool enable_rt_engine = Dune::Dorie::RaviartThomasFluxReconstructionHelper<dim,flux_order,BaseTraits::GridGeometryType>::enable;
225
  using GFFluxReconstruction = std::conditional_t<enable_rt_engine,Dune::Dorie::RaviartThomasFluxReconstruction<GV,RF,flux_order,BaseTraits::GridGeometryType>,void>;
226 227 228 229
  /// Peclet
  using GFPeclet = Dune::Dorie::PecletAdapter<BaseTraits,SoluteParameters,GFWaterFlux,GFWaterContent>;
  /// Effective Hydrodynamic Dipspersion
  using GFEffectiveHydrodynamicDispersion = Dune::Dorie::EffectiveHydrodynamicDispersionAdapter<BaseTraits,SoluteParameters,GFWaterFlux,GFWaterContent>;
230

231 232 233
  // -- Utility Class Definitions -- //
  /// VTK Output writer base class
  using Writer = Dune::Dorie::GridFunctionVTKSequenceWriter<GV>;
234 235 236 237 238 239 240 241

  /// Policy when checking for negative values in the solution
  enum class CheckNegativePolicy
  {
    Pass, //!< Always pass the check
    Warn, //!< Pass the check but issue warnings for negative values
    Retry //!< Reject the solution for negative values
  };
242 243 244 245
};


/*-------------------------------------------------------------------------*//**
246 247 248
 * @brief      Class for Transport model. 
 * @details    For given ModelTransportTraits, this class manages the 
 *             model for the transport equation.
249 250 251
 * @author     Santiago Ospina De Los Ríos
 * @date       2018
 * @ingroup    TransportModel
252
 * @ingroup    Collective
253
 *
254 255 256 257
 * @todo  Implement source term.
 * @todo  Implement more complex initial conditions.
 * @todo  Implement adaptivity.
 * @todo  Implement DG local operator.
258
 * 
259
 * @tparam     ModelTransportTraits  Traits containing the type definitions
260 261
 *                                        which this class will use.
 */
262 263
template<typename ModelTransportTraits>
class ModelTransport : public ModelBase
264 265
{
public:
266
  using Traits = ModelTransportTraits;
267 268 269 270

protected:

  static constexpr int dim = Traits::dim;
271
  static constexpr int order = Traits::order;
272 273 274
  static constexpr int flux_order = Traits::flux_order;
  static constexpr bool enable_rt_engine = Traits::enable_rt_engine;

275 276 277 278 279 280 281 282
  using TimeField = typename Traits::TimeField;
  using RF = typename Traits::RF;
  using Grid = typename Traits::Grid;
  using GV = typename Traits::GV;

  /// (Inestationary) Grid Function Type for water flux
  using GFWaterFlux = typename Traits::GFWaterFlux;
  /// (Inestationary) Grid Function Type for water content
283
  using GFWaterContent = typename Traits::GFWaterContent;
284 285 286 287 288 289 290 291

  /// GFS Helper
  using GFSHelper = typename Traits::GFSHelper;
  /// Problem GFS
  using GFS = typename Traits::GFS;
  /// GFS Constraints Container
  using CC = typename Traits::CC;

292 293 294 295 296 297 298
  /// LSGFS Helper
  using LSGFSHelper = typename Traits::LSGFSHelper;
  /// Linear solver GFS
  using LSGFS = typename Traits::LSGFS;
  /// LSGFS Constraints Container
  using LSCC = typename Traits::LSCC;

299
  // -- D/*ORiE Class Definitions -- //
300
  /// Wrapper for grid and volume & boundary mappings
301
  using GridCreator = typename Traits::GridCreator;
302 303
  /// Class defining the solute parameterization
  using SoluteParameters = typename Traits::SoluteParameters;
304 305 306 307
  /// Class defining boundary conditions
  using SoluteBoundary = typename Traits::SoluteBoundary;
  /// Class defining the initial condition
  using SoluteInitial = typename Traits::SoluteInitial;
308 309
  /// Class defining the initial condition factory
  using SoluteInitialFactory = typename Traits::SoluteInitialFactory;
310 311 312 313 314
  /// Local spatial operator
  using SLOP = typename Traits::SLOP;
  /// Local temporal operator
  using TLOP = typename Traits::TLOP;
  /// Time controller
315
  using TimeStepController = typename Traits::TimeStepController;
316 317 318 319 320 321 322 323 324

  // -- DUNE Class Defintions -- //
  /// Matrix backend type
  using MBE = typename Traits::MBE;
  /// Grid operator for spatial LOP
  using GO0 = typename Traits::GO0;
  /// Grid operator for temporal LOP
  using GO1 = typename Traits::GO1;
  /// One step grid operator for merging both GOs
325 326
  using ImplicitIGO = typename Traits::ImplicitIGO;
  using ExplicitIGO = typename Traits::ExplicitIGO;
327 328 329
  /// Solution vector type
  using U = typename Traits::U;
  /// Linear solver type
330 331
  using ImplicitLS = typename Traits::ImplicitLS;
  using ExplicitLS = typename Traits::ExplicitLS;
332
  /// Methods computing the time step
333 334 335
  using ImplicitSLPS = typename Traits::ImplicitSLPS;
  using ImplicitOSM = typename Traits::ImplicitOSM;
  using ExplicitOSM = typename Traits::ExplicitOSM;
336

337
  /// Model state provider
338
  using State = typename Traits::State;
339
  using ConstState = typename Traits::ConstState;
340

341 342
  // -- Grid Functions of state variables -- //
  /// Solute concentration
343
  using GFSolute = typename Traits::GFSolute;
344
  /// Total solute
345
  using GFTotalSolute = typename Traits::GFTotalSolute;
346 347
  /// Water Flux reconstruction
  using GFFluxReconstruction = typename Traits::GFFluxReconstruction;
348 349 350 351
  /// Peclet
  using GFPeclet = typename Traits::GFPeclet;
  /// Effective Hydrodynamic Dipspersion
  using GFEffectiveHydrodynamicDispersion = typename Traits::GFEffectiveHydrodynamicDispersion;
352

353 354 355 356 357
  // // -- Utility Class Definitions -- //
  /// VTK Output writer base class
  using Writer = Dune::Dorie::GridFunctionVTKSequenceWriter<GV>;

  /// Time stepping scheme
358
  using TimeSteppingParameter = Dune::PDELab::TimeSteppingParameterInterface<RF>;
359

360 361 362
  /// Policy for negative values in solution
  using CheckNegativePolicy = typename Traits::CheckNegativePolicy;

363 364

  // -- Member variables -- //
365 366
  Dune::MPIHelper& helper;
  Dune::ParameterTree& inifile;
367

368 369 370
  /// Time step controller of this instance
  TimeStepController time_step_ctrl;

371 372 373
  const Dune::VTK::OutputType output_type;

  std::shared_ptr<Grid> grid;
374 375 376 377
  GV gv;
  std::shared_ptr<GFS> gfs;
  std::unique_ptr<CC> cc;

378
  std::shared_ptr<SoluteParameters> sparam;
379
  std::shared_ptr<SoluteBoundary> sboundary;
380
  std::shared_ptr<SoluteInitial> sinitial;
381

382 383 384
  std::shared_ptr<LSGFS> lsgfs;
  std::shared_ptr<LSCC> lscc;

385 386 387 388
  std::unique_ptr<SLOP> slop;
  MBE mbe_slop;
  std::unique_ptr<TLOP> tlop;
  MBE mbe_tlop;
389
  std::shared_ptr<GO0> go0;
390
  std::unique_ptr<GO1> go1;
391 392 393 394 395 396 397
  std::unique_ptr<ImplicitIGO> implicit_igo;
  std::unique_ptr<ExplicitIGO> explicit_igo;
  std::unique_ptr<ImplicitLS>  implicit_ls;
  std::unique_ptr<ExplicitLS>  explicit_ls;
  std::unique_ptr<ImplicitSLPS> implicit_slps;
  std::unique_ptr<ImplicitOSM> implicit_osm;
  std::unique_ptr<ExplicitOSM> explicit_osm;
398
  std::unique_ptr<TimeSteppingParameter> ts_param;
399 400 401

  std::shared_ptr<U> u;

402 403
  std::shared_ptr<GFWaterFlux>            water_flux_gf;
  std::shared_ptr<GFWaterContent>         water_content_gf;
404 405 406

  std::unique_ptr<Writer> vtkwriter;

407 408
  mutable std::shared_ptr<GFFluxReconstruction> cache_solute_flux_gf_rt;

409 410
  unsigned int time_steps;

411 412
  const RF _courant_number;

413
  const bool enable_fluxrc;
414

415 416 417
  /// Stored policy for negative values in solution
  CheckNegativePolicy check_negative_policy;

418 419
public:

420
  /*-----------------------------------------------------------------------*//**
421
   * @brief      Construct model.
422
   *
423 424 425
   * @param      _inifile           Dune parameter file parser.
   * @param      _grid_creator      Mapper for grid and volume/boundary data.
   * @param      _helper            Dune MPI instance controller.
426
   */
427
  ModelTransport (
428 429
    Dune::ParameterTree&            _inifile,
    const GridCreator&              _grid_creator,
430
    Dune::MPIHelper&                _helper
431 432
  );

433
  ~ModelTransport() = default;
434

435 436
  /*------------------------------------------------------------------------*//**
   * @brief      Compute a time step. Models the time step requirement of
437
   *             ModelBase
438
   *
439
   * @throws     Dune::Exception  Fatal error occurs during computation.
440 441 442
   */
  void step () override;

443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458
  /*------------------------------------------------------------------------*//**
   * @brief      Mark the grid in order to improve the current model.
   */
  void mark_grid() override;

  /*------------------------------------------------------------------------*//**
   * @brief      Adapt the grid and the solutions of the transport simulation
   * @warning    This method does not adapt solutions of other simulations!
   */
  void adapt_grid() override;

  /*------------------------------------------------------------------------*//**
   * @brief      Setup operators to fit the new grid.
   */
  void post_adapt_grid() override;

459 460 461 462 463
  /*-----------------------------------------------------------------------*//**
   * @brief      Method that provides the begin time of the model.
   *
   * @return     Begin time of the model.
   */
464
  double begin_time() const final {return time_step_ctrl.get_time_begin();}
465 466 467 468 469 470

  /*-----------------------------------------------------------------------*//**
   * @brief      Method that provides the end time of the model.
   *
   * @return     End time of the model.
   */
471
  double end_time() const final {return time_step_ctrl.get_time_end();}
472

473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492
  /*------------------------------------------------------------------------*//**
   * @brief      Provides an object with all the information that describes
   *             the current state of the solved variable.
   *
   * @return     Current state of the model.
   */
  State current_state() {
    return State{gfs,u,current_time()};
  }

  /*------------------------------------------------------------------------*//**
   * @brief      Provides an object with all the information that describes
   *             the current state of the solved variable.
   *
   * @return     (Const) Current state of the model.
   */
  ConstState current_state() const {
    return ConstState{gfs,u,current_time()};
  }

493 494 495 496 497 498
  /*-----------------------------------------------------------------------*//**
   * @brief      Suggest a time step to the model.
   *
   * @param[in]  dt    Suggestion for the internal time step of the model. The
   *                   new internal time step shall not be bigger than dt.
   */
499 500 501 502
  void suggest_timestep(double dt) final
  {
    time_step_ctrl.suggest_time_step(dt);
  }
503

504 505 506 507 508 509
  /*------------------------------------------------------------------------*//**
   * @brief      Write the current data using the VTKWriter. VTKWriters is 
   *             always cleared after used.
   */
  virtual void write_data () const override;

510
  /*-----------------------------------------------------------------------*//**
511
   * @brief      Sets an instationary function for the water flux.
512
   *
513
   * @param[in]  _water_flux_gf  The grid function water flux.
514
   */
515
  void set_water_flux(std::shared_ptr<GFWaterFlux> _water_flux_gf)
516
  {
517
    if (_water_flux_gf)
518
      water_flux_gf = _water_flux_gf;
519
    else
520
      DUNE_THROW(Dune::InvalidStateException,
521
        "Pointer to water_flux_gf is invalid!");
522

523
    slop->set_water_flux(_water_flux_gf);
524 525
  }

526
  /*-----------------------------------------------------------------------*//**
527
   * @brief      Sets an instationary function for the water content.
528
   *
529
   * @param[in]  _water_content_gf  The grid function water content.
530
   */
531
  void set_water_content(std::shared_ptr<GFWaterContent> _water_content_gf)
532
  {
533
    if (_water_content_gf)
534
      water_content_gf = _water_content_gf;
535
    else
536
      DUNE_THROW(Dune::InvalidStateException,
537
        "Pointer to water_content_gf is invalid!");
538

539 540
    slop->set_water_content(_water_content_gf);
    tlop->set_water_content(_water_content_gf);
541 542
  }

543 544
  /*------------------------------------------------------------------------*//**
   * @brief      Generates a solute grid function from a given state
545
   *
546
   * @return     Pointer to a solute grid function
547
   */
548
  std::shared_ptr<const GFSolute> get_solute(ConstState state) const
549
  {
550
    return std::make_shared<GFSolute>(state.grid_function_space,state.coefficients);
551 552 553
  }

  /*------------------------------------------------------------------------*//**
554 555 556
   * @brief      Gets the solute grid function
   *
   * @return     Pointer to a solute grid function
557
   */
558
  std::shared_ptr<const GFSolute> get_solute() const
559 560 561
  {
    return get_solute(current_state());
  }
562

563 564 565 566 567 568 569 570 571
  /*------------------------------------------------------------------------*//**
   * @brief      Generates a total solute grid function from a given state
   *
   * @return     Pointer to a total solute grid function
   */
  std::shared_ptr<const GFTotalSolute> get_total_solute(ConstState state) const
  {
    auto solute_gf = get_solute(state);

572 573 574 575 576
    auto _water_content_gf 
      = std::make_shared<GFWaterContent>(*water_content_gf);
    _water_content_gf->setTime(state.time);
    _water_content_gf->shrink_to_time();
    
577
    return std::make_shared<GFTotalSolute>(solute_gf, _water_content_gf);
578 579 580 581 582 583 584 585 586 587 588 589
  }

  /*------------------------------------------------------------------------*//**
   * @brief      Gets the total solute grid function
   *
   * @return     Pointer to a total solute grid function
   */
  std::shared_ptr<const GFTotalSolute> get_total_solute() const
  {
    return get_total_solute(current_state());
  }

590
  /*-----------------------------------------------------------------------*//**
591 592 593 594
   * @brief      Generates a (reconstructed) solute flux grid function from a
   *             given state
   *
   * @param[in]  state  The state
595
   *
596
   * @return     Pointer to a solute grid function
597
   */
598 599 600
  template<bool enabled=enable_rt_engine>
	std::enable_if_t<enabled, std::shared_ptr<const GFFluxReconstruction>>
  get_solute_flux_reconstructed(ConstState state) const
601
  {
602
    std::shared_ptr<GFFluxReconstruction> gf_ptr;
603 604 605

    auto& cache = cache_solute_flux_gf_rt;

606 607 608
    if (state.grid_function_space != gfs
        or state.coefficients != u
        or state.time != current_time())
609
    {
610
      // if state is different to current state, create flux from zero
611

612 613
      gf_ptr = std::make_unique<GFFluxReconstruction>(
                      this->_log,gv,inifile.sub("fluxReconstruction"));
614

615
      slop->setTime(state.time);
616

617 618 619 620
      // update it with the state
      gf_ptr->update(*(state.coefficients),
                      *(state.grid_function_space),
                      *slop);
621

622 623 624
      slop->setTime(current_time());
    } else if (not cache) {
      // if state is equal to current state, use cache.
625

626 627
      cache = std::make_unique<GFFluxReconstruction>(
                      this->_log,gv,inifile.sub("fluxReconstruction"));
628

629
      slop->setTime(state.time);
630

631 632
      // update it with current state
      cache->update(*u,*gfs,*slop);
633

634
      slop->setTime(current_time());
635

636
      gf_ptr = cache;
637
    } else {
638
      gf_ptr = cache;
639
    }
640

641
    return gf_ptr;
642 643 644
  }

  /*------------------------------------------------------------------------*//**
645 646 647
   * @brief      Gets the (reconstructed) solute flux grid function
   *
   * @return     Pointer to a solute grid function
648
   */
649 650 651
  template<bool enabled=enable_rt_engine>
	std::enable_if_t<enabled, std::shared_ptr<const GFFluxReconstruction>>
  get_solute_flux_reconstructed() const
652 653 654
  {
    return get_solute_flux_reconstructed(current_state());
  }
655 656 657 658 659 660 661
protected:

  /*------------------------------------------------------------------------*//**
   * @brief      Create the Grid Operators, Solvers, and Time Step Operator.
   *             This function must be called after changes to grid or GFS!
   */
  void operator_setup ();
662

663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679
  /// Check the solution of a time step
  /** This function currently performs the following checks:
   *
   *  * No negative values in the resulting solution \p u_after
   *
   *  The method returns a pair of `bool` and `string`. The boolean indicates
   *  if the check was successful, and the string contains a check error
   *  message if it was not.
   *
   *  \param u_before Solution at the start of the time step
   *  \param u_after Computed solution at the end of the time step
   *  \return Pair of `bool` (check passed) and `string` (check fail
   *          message, if any)
   */
  std::pair<bool, std::string> check_solution (const U& u_before,
                                               const U& u_after) const;

680 681 682 683 684 685
private:

#ifdef GTEST
  FRIEND_TEST(ModelCoupling, TimeSteps);
#endif

686 687 688 689 690
};

} // namespace Dorie
} // namespace Dune

691
#endif // DUNE_DORIE_MODEL_TRANSPORT_HH