richards.hh 21.4 KB
Newer Older
1 2
#ifndef DUNE_DORIE_MODEL_RICHARDS_HH
#define DUNE_DORIE_MODEL_RICHARDS_HH
3

4 5 6 7 8 9 10 11 12 13 14 15
#include <memory>
#include <string>
#include <iostream>

#include <dune/common/parametertree.hh>
#include <dune/common/parallel/mpihelper.hh>

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

16 17 18
#include <dune/dorie/common/util.hh>
#include <dune/dorie/common/interpolator.hh>
#include <dune/dorie/common/time_controller.hh>
19
#include <dune/dorie/common/grid_creator.hh>
20
#include <dune/dorie/common/grid_function_writer.hh>
21
#include <dune/dorie/common/flux_reconstruction/raviart_thomas.hh>
22

23
#include <dune/dorie/model/base.hh>
24
#include <dune/dorie/model/richards/initial_condition.hh>
25
#include <dune/dorie/model/richards/adaptivity/handler.hh>
26
#include <dune/dorie/model/richards/parameterization/interface.hh>
27 28 29 30
#include <dune/dorie/model/richards/adapters/water_content.hh>
#include <dune/dorie/model/richards/adapters/saturation.hh>
#include <dune/dorie/model/richards/adapters/water_flux.hh>
#include <dune/dorie/model/richards/adapters/conductivity.hh>
31
#include <dune/dorie/model/richards/flow_parameters.hh>
32 33
#include <dune/dorie/model/richards/flow_boundary.hh>
#include <dune/dorie/model/richards/flow_source.hh>
34
#include <dune/dorie/model/richards/local_operator_DG.hh>
35
#include <dune/dorie/model/richards/local_operator_FV.hh>
36
#include <dune/dorie/model/richards/local_operator_cfg.hh>
37 38 39 40 41

namespace Dune{
namespace Dorie{

/*-------------------------------------------------------------------------*//**
42
 * @brief      Traits for the class ModelRichards. In particular, this
43
 *             class extends BaseTraits to be used for the DG implementation of
44
 *             the Richards model
45
 *
46
 * @ingroup    ModelRichards
47
 *
48
 * @tparam     BaseTraits  Traits defining domain and range field properties of
49
 *                         the model.
50
 * @tparam     k           Order of the polynomial degree used for the basis
51 52
 *                         functions in the DG method
 */
53
template<class BaseTraits, int k>
54
struct ModelRichardsTraits : public BaseTraits
55 56
{
	static constexpr int dim = BaseTraits::dim;
57
	static constexpr int order = k;
58
	static constexpr int flux_order = (order==0) ? 0 : (order-1);
59
	using RF = typename BaseTraits::RF;
60
	using Scalar = typename BaseTraits::Scalar;
61 62 63
	using Grid = typename BaseTraits::Grid;
	using GV = typename BaseTraits::GV;

64 65 66 67 68 69
  /// GFS Helper
  using GFSHelper = 
    std::conditional_t<order==0,
      Dune::Dorie::LinearSolverGridFunctionSpaceHelper<GV,RF,BaseTraits::GridGeometryType>,
      Dune::Dorie::GridFunctionSpaceHelper<GV,RF,order,BaseTraits::GridGeometryType>>;

70 71 72 73 74 75 76 77 78 79 80 81 82 83
	/// Problem GFS
	using GFS = typename GFSHelper::Type;
	/// GFS Constraints Container
	using CC = typename GFSHelper::CC;


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

	// -- DORiE Class Definitions -- //
84
	/// Wrapper for grid and volume & boundary mappings
85
	using GridCreator = Dune::Dorie::GridCreator<typename BaseTraits::Grid>;
86
	/// Class defining the soil parameterization
87
	using FlowParameters = Dune::Dorie::FlowParameters<BaseTraits>;
88 89 90
	/// Class defining boundary conditions
	using FlowBoundary = Dune::Dorie::FlowBoundary<BaseTraits>;
	/// Class defining source terms
91
	using FlowSource = Dune::Dorie::FlowSource<BaseTraits>;
92
	/// Class defining the initial condition
93
	using FlowInitial = Dune::Dorie::InitialCondition<BaseTraits>;
94
	/// Class defining the initial condition factory
95
	using FlowInitialFactory = Dune::Dorie::RichardsInitialConditionFactory<BaseTraits>;
96
	/// Local spatial operator
97 98 99 100 101
  using SLOP = 
    std::conditional_t<order==0,
      Dune::Dorie::Operator::RichardsFVSpatialOperator<FlowParameters,FlowBoundary>,
      Dune::Dorie::Operator::RichardsDGSpatialOperator<BaseTraits,FlowParameters,FlowBoundary,FlowSource,typename GFSHelper::FEM,false>>;

102
	/// Local temporal operator
103 104 105 106 107
  using TLOP = 
    std::conditional_t<order==0,
      Dune::Dorie::Operator::RichardsFVTemporalOperator<FlowParameters>,
      Dune::Dorie::Operator::RichardsDGTemporalOperator<BaseTraits,FlowParameters,typename GFSHelper::FEM,false>>;

108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131
	/// Time controller
	using CalculationController = Dune::Dorie::CalculationController<RF,FlowBoundary>;

	// -- 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>;
	/// One step grid operator for merging both GOs
	using IGO = Dune::PDELab::OneStepGridOperator<GO0,GO1>;
	/// Solution vector type
	using U = typename IGO::Traits::Domain;
	/// Linear solver types
	using LS = Dune::PDELab::ISTLBackend_OVLP_AMG_4_DG<IGO,CC,LSGFS,LSCC,
		Dune::PDELab::CG2DGProlongation,Dune::SeqSSOR,Dune::BiCGSTABSolver>;
	/// Non-linear solver types
	using PDESOLVER = Dune::PDELab::Newton<IGO,LS,U>;
	/// Time stepping scheme
	using TimeStepScheme = Dune::PDELab::Alexander2Parameter<RF>;
	/// Methods computing the time step
	using OSM = Dune::PDELab::OneStepMethod<RF,IGO,PDESOLVER,U,U>;

132
	/// Model state provider
133 134 135 136 137 138 139 140 141 142 143
	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;
	};

144
	/// Model state provider
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
	struct State 
	{
		using GridFunctionSpace = GFS;
		using Coefficients = U;
		using TimeField = typename BaseTraits::TimeField;

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

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

161
	// -- Grid Functions of state variables -- //
162
	/// Matric head
163
	using GFMatricHead = Dune::PDELab::DiscreteGridFunction<GFS,U>;
164
	/// Water Flux
165
	using GFWaterFlux = Dune::Dorie::WaterFluxAdapter<GFS,U,FlowParameters>;
166
	/// Conductivity
167
	using GFConductivity = Dune::Dorie::ConductivityAdapter<BaseTraits,FlowParameters>;
168
	/// Water content
169
	using GFWaterContent = Dune::Dorie::WaterContentAdapter<BaseTraits,FlowParameters,GFMatricHead>;
170
	/// Saturation
171
	using GFSaturation = Dune::Dorie::SaturationAdapter<BaseTraits,FlowParameters,GFMatricHead>;
172
	/// Water Flux reconstruction
173
	static constexpr bool enable_rt_engine = Dune::Dorie::RaviartThomasFluxReconstructionHelper<dim,flux_order,BaseTraits::GridGeometryType>::enable;
174
	using GFFluxReconstruction = std::conditional_t<enable_rt_engine,Dune::Dorie::RaviartThomasFluxReconstruction<GV,RF,flux_order,BaseTraits::GridGeometryType>,void>;
175

176
	// -- Utility Class Definitions -- //
177
	/// Custom VTK output writer
Santiago Ospina's avatar
Santiago Ospina committed
178
	using Writer = Dune::Dorie::GridFunctionVTKSequenceWriter<GV>;
179
	/// Grid Adaptivity Handler base class
180
	using AdaptivityHandler = Dune::Dorie::AdaptivityHandlerBase<BaseTraits,GFS,FlowParameters,FlowBoundary,U,order>;
181
	/// Factory for creating the AdaptivityHandler
182
	using AdaptivityHandlerFactory = Dune::Dorie::AdaptivityHandlerFactory<BaseTraits,GFS,FlowParameters,FlowBoundary,U,order>;
183 184 185
};

/*-------------------------------------------------------------------------*//**
186
 * @brief      Class for Richards model. This class manages the model
187 188 189
 *             for the richards equation. It can perform time-steps and print
 *             the solution
 *
190
 * @ingroup    ModelRichards
191
 *
192 193
 * @tparam     ModelRichardsTraits  Traits containing the type definitions
 *                                       for Richards Model
194
 */
195 196
template<typename ModelRichardsTraits>
class ModelRichards : public ModelBase
197 198
{
public:
199
	using Traits = ModelRichardsTraits;
200 201 202 203

private:

	static constexpr int dim = Traits::dim;
204 205
	static constexpr int order = Traits::order;
	static constexpr int flux_order = Traits::flux_order;
206
	static constexpr bool enable_rt_engine = Traits::enable_rt_engine;
207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
	using RF = typename Traits::RF;
	using Grid = typename Traits::Grid;
	using GV = typename Traits::GV;

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


	/// LSGFS Helper
	using LSGFSHelper = typename Traits::LSGFSHelper;
	/// Linear solver GFS
	using LSGFS = typename Traits::LSGFS;
	/// LSGFS Constraints Container
	using LSCC = typename Traits::LSCC;

	// -- DORiE Class Definitions -- //
227
	/// Wrapper for grid and volume & boundary mappings
228
	using GridCreator = typename Traits::GridCreator;
229
	/// Class defining the soil parameterization
230
	using FlowParameters = typename Traits::FlowParameters;
231 232 233 234 235 236
	/// Class defining boundary conditions
	using FlowBoundary = typename Traits::FlowBoundary;
	/// Class defining source terms
	using FlowSource = typename Traits::FlowSource;
	/// Class defining the initial condition
	using FlowInitial = typename Traits::FlowInitial;
237 238
	/// Class defining the initial condition factory
	using FlowInitialFactory = typename Traits::FlowInitialFactory;
239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266
	/// Local spatial operator
	using SLOP = typename Traits::SLOP;
	/// Local temporal operator
	using TLOP = typename Traits::TLOP;
	/// Time controller
	using CalculationController = typename Traits::CalculationController;

	// -- 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
	using IGO = typename Traits::IGO;
	/// Solution vector type
	using U = typename Traits::U;
	/// Linear solver types
	using LS = typename Traits::LS;

	/// Non-linear solver types
	using PDESOLVER = typename Traits::PDESOLVER;
	/// Time stepping scheme
	using TimeStepScheme = typename Traits::TimeStepScheme;
	/// Methods computing the time step
	using OSM = typename Traits::OSM;

267
	/// Model state provider
268 269 270
	using State = typename Traits::State;
	using ConstState = typename Traits::ConstState;

271 272 273 274 275 276 277 278 279 280 281
	// -- Grid Functions of state valriables -- //
	/// Matric head
	using GFMatricHead = typename Traits::GFMatricHead;
	/// Water Flux
	using GFWaterFlux = typename Traits::GFWaterFlux;
	/// Conductivity
	using GFConductivity = typename Traits::GFConductivity;
	/// Water content
	using GFWaterContent = typename Traits::GFWaterContent;
	/// Saturation
	using GFSaturation = typename Traits::GFSaturation;
282
	/// Water Flux reconstruction
283
	using GFFluxReconstruction = typename Traits::GFFluxReconstruction;
284 285


286 287
	// -- Utility Class Definitions -- //
	/// VTK Output writer base class
288
	using Writer = typename Traits::Writer;
289 290 291 292 293 294 295
	/// Grid Adaptivity Handler base class
	using AdaptivityHandler = typename Traits::AdaptivityHandler;
	/// Factory for creating the AdaptivityHandler
	using AdaptivityHandlerFactory = typename Traits::AdaptivityHandlerFactory;

protected:

296
	// -- Member variables -- //
297
	Dune::MPIHelper& helper;
298
	Dune::ParameterTree inifile;
299

300 301
	const Dune::VTK::OutputType output_type;

302
	std::shared_ptr<GFS> gfs;
303
	std::shared_ptr<Grid> grid;
304
	GV gv;
305 306
	std::unique_ptr<CC> cc;

307 308
	std::shared_ptr<LSGFS> lsgfs;
	std::shared_ptr<LSCC> lscc;
309

310
	std::shared_ptr<FlowParameters> fparam;
311 312
	std::shared_ptr<const FlowBoundary> fboundary;
	std::shared_ptr<const FlowSource> fsource;
313 314 315 316 317 318 319 320 321 322 323 324 325 326 327
	std::unique_ptr<FlowInitial> finitial;
	std::unique_ptr<CalculationController> controller;

	std::unique_ptr<SLOP> slop;
	MBE mbe_slop;
	std::unique_ptr<TLOP> tlop;
	MBE mbe_tlop;
	std::unique_ptr<GO0> go0;
	std::unique_ptr<GO1> go1;
	std::unique_ptr<IGO> igo;
	std::unique_ptr<LS> ls;
	std::unique_ptr<PDESOLVER> pdesolver;
	std::unique_ptr<OSM> osm;
	TimeStepScheme alex2;

328
	std::shared_ptr<U> u;
329 330

	std::unique_ptr<Writer> vtkwriter;
331 332
	std::unique_ptr<AdaptivityHandler> adaptivity;

333 334
	mutable std::shared_ptr<GFFluxReconstruction> cache_water_flux_gf_rt;

335 336
	RF time_before;
	RF dt_before;
337

338
	const bool enable_fluxrc;
339 340
public:

341
	/*------------------------------------------------------------------------*//**
342
	 * @brief      Construct Richards model.
343
	 *
344 345 346
	 * @param      _inifile      Dune parameter file parser
	 * @param      _grid_mapper  Mapper for grid and volume/boundary data
	 * @param      _helper       Dune MPI instance controller
347
	 */
348 349
	ModelRichards (
		const Dune::ParameterTree& _inifile,
350
		const GridCreator& _grid_creator,
351 352
		Dune::MPIHelper& _helper
	);
353

354
	~ModelRichards() = default;
355

356
	/*------------------------------------------------------------------------*//**
357
	 * @brief      Compute a time step. Models the time step requirement of
358
	 *             ModelBase
359 360
	 *
	 * @throws     Dune::Exception  Fatal error occurs during computation
361
	 */
362
	void step () override;
363

364
	/*------------------------------------------------------------------------*//**
365 366 367
	 * @brief      Mark the grid in order to improve the current model.
	 */
	virtual void mark_grid() override;
368

369
	/*------------------------------------------------------------------------*//**
370
	 * @brief      Adapt the grid together it every dependency of the
371 372 373 374
	 *             grid (e.g. solution vector and grid function spaces).
	 */
	virtual void adapt_grid() override;

375
	/*------------------------------------------------------------------------*//**
376 377 378 379
	 * @brief      Setup operators to fit the new grid.
	 */
	virtual void post_adapt_grid() override;

380 381
	/*------------------------------------------------------------------------*//**
	 * @brief      Method that provides the begin time of the model.
382
	 *
383 384
	 * @return     Begin time of the model.
	 */
385
	double begin_time() const final {return controller->getBeginTime();}
386

387 388
	/*------------------------------------------------------------------------*//**
	 * @brief      Method that provides the end time of the model.
389
	 *
390
	 * @return     End time of the model.
391
	 */
392
	double end_time() const final {return controller->getEndTime();}
393

394 395 396 397 398
	/*------------------------------------------------------------------------*//**
	 * @brief      Method that provides the current time of the model.
	 *
	 * @return     Current time of the model.
	 */
399
	double current_time() const final {return controller->getTime();}
400

401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420
	/*------------------------------------------------------------------------*//**
	 * @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()};
	}

421
	/*------------------------------------------------------------------------*//**
422 423
	 * @brief      Suggest a time step to the model.
	 *
424
	 * @param[in]  dt    Suggestion for the internal time step of the model. The
425
	 *                   new internal time step shall not be bigger than dt.
426
	 */
427
	void suggest_timestep(double dt) final {controller->suggest_timestep(dt);}
428 429 430 431 432

	/*------------------------------------------------------------------------*//**
	 * @brief      Write the data using the VTKWriter. VTKWriters is always
	 *             cleared after used.
	 */
433 434
	void virtual write_data() const override;

435
	/*------------------------------------------------------------------------*//**
436 437 438 439
	 * @brief      Generates a matric head grid function from a given state
	 *
	 * @return     Pointer to a matric head grid function
	 */
440
	std::shared_ptr<const GFMatricHead> get_matric_head(ConstState state) const
441
	{
442
		return std::make_shared<GFMatricHead>(state.grid_function_space,state.coefficients);
443
	}
444

445 446
	/*------------------------------------------------------------------------*//**
	 * @brief      Gets the matric head grid function
447
	 *
448
	 * @return     Pointer to a matric head grid function
449
	 */
450
	std::shared_ptr<const GFMatricHead> get_matric_head() const
451
	{
452
		return get_matric_head(current_state());
453
	}
454

455
	/*------------------------------------------------------------------------*//**
456
	 * @brief      Generates a conductivity grid function from a given state
457
	 *
458
	 * @return     Pointer to a conductivity grid function
459
	 */
460
	std::shared_ptr<const GFConductivity> get_conductivity(ConstState state) const
461
	{
462 463 464 465 466
		// hard copy of parameters
		auto _fparam = std::make_shared<FlowParameters>(*fparam);

		// create conductivity adapter
		return std::make_shared<GFConductivity>(gv, _fparam);
467 468 469 470 471 472 473
	}

	/*------------------------------------------------------------------------*//**
	 * @brief      Gets the conductivity grid function
	 *
	 * @return     Pointer to a conductivity grid function
	 */
474
	std::shared_ptr<const GFConductivity> get_conductivity() const
475 476 477
	{
		return get_conductivity(current_state());
	}
478

479 480 481 482 483
	/*------------------------------------------------------------------------*//**
	 * @brief      Generates a water content grid function from a given state
	 *
	 * @return     Pointer to a water content grid function
	 */
484
	std::shared_ptr<const GFWaterContent> get_water_content(ConstState state) const
485
	{
486 487 488 489 490
		// hard copy of parameters
		auto _fparam = std::make_shared<FlowParameters>(*fparam);

		// create water content adapter
		return std::make_shared<GFWaterContent>(get_matric_head(state), gv, _fparam);
491 492 493 494 495 496 497
	}

	/*------------------------------------------------------------------------*//**
	 * @brief      Gets the water content grid function
	 *
	 * @return     Pointer to a water content grid function
	 */
498
	std::shared_ptr<const GFWaterContent> get_water_content() const
499 500 501 502 503 504 505 506 507
	{
		return get_water_content(current_state());
	}

	/*------------------------------------------------------------------------*//**
	 * @brief      Generates a saturation grid function from a given state
	 *
	 * @return     Pointer to a saturation grid function
	 */
508
	std::shared_ptr<const GFSaturation> get_saturation(ConstState state) const
509
	{
510 511 512 513 514
		// hard copy of parameters
		auto _fparam = std::make_shared<FlowParameters>(*fparam);

		// create saturation adapter
		return std::make_shared<GFSaturation>(get_matric_head(state), gv, _fparam);
515 516 517 518 519 520 521
	}

	/*------------------------------------------------------------------------*//**
	 * @brief      Gets the saturation grid function
	 *
	 * @return     Pointer to a saturation grid function
	 */
522
	std::shared_ptr<const GFSaturation> get_saturation() const
523 524 525 526 527 528 529 530 531
	{
		return get_saturation(current_state());
	}

	/*------------------------------------------------------------------------*//**
	 * @brief      Generates a water flux grid function from a given state
	 *
	 * @return     Pointer to a water flux grid function
	 */
532
	std::shared_ptr<const GFWaterFlux> get_water_flux(ConstState state) const
533
	{
534 535 536 537 538
		// hard copy of parameters
		auto _fparam = std::make_shared<FlowParameters>(*fparam);

		// create water flux adapter
		return std::make_shared<GFWaterFlux>(state.grid_function_space, state.coefficients, _fparam);
539 540 541 542 543 544 545
	}

	/*------------------------------------------------------------------------*//**
	 * @brief      Gets the water flux grid function
	 *
	 * @return     Pointer to a water flux grid function
	 */
546
	std::shared_ptr<const GFWaterFlux> get_water_flux() const
547 548 549 550 551 552 553 554 555 556 557 558
	{
		return get_water_flux(current_state());
	}

	/*------------------------------------------------------------------------*//**
	 * @brief      Generates a (reconstructed) water flux grid function from a
	 *             given state
	 *
	 * @param[in]  state  The state
	 *
	 * @return     Pointer to a (reconstructed) water flux grid function
	 */
559 560 561
	template<bool enabled=enable_rt_engine>
	std::enable_if_t<enabled, std::shared_ptr<const GFFluxReconstruction>>
	get_water_flux_reconstructed(ConstState state) const
562
	{
563 564 565 566
		std::shared_ptr<GFFluxReconstruction> gf_ptr;

		auto& cache = cache_water_flux_gf_rt;

567 568 569
		if (state.grid_function_space != gfs
			or state.coefficients != u
			or state.time != current_time())
570
		{
571
			// if state is different to current state, create flux from zero
572

573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595
			gf_ptr = std::make_unique<GFFluxReconstruction>(
									this->_log,gv,inifile.sub("fluxReconstruction"));

			slop->setTime(state.time);

			// update it with the state
			gf_ptr->update(*(state.coefficients),
										*(state.grid_function_space),
										*slop);

			slop->setTime(current_time());
		}

		// if state is equal to current state, use cache.
		else if (not cache) {

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

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

			gf_ptr = cache;
596
		}
597 598 599 600 601

		else {
			gf_ptr = cache;
		}

602
		return gf_ptr;
603 604 605 606 607 608 609
	}

	/*------------------------------------------------------------------------*//**
	 * @brief      Gets the (reconstructed) water flux grid function
	 *
	 * @return     Pointer to a (reconstructed) water flux grid function
	 */
610 611 612
	template<bool enabled=enable_rt_engine>
	std::enable_if_t<enabled, std::shared_ptr<const GFFluxReconstruction>>
	get_water_flux_reconstructed() const
613 614
	{
		return get_water_flux_reconstructed(current_state());
615
	}
616

617 618 619
protected:

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

625 626 627 628 629
};

} // namespace Dorie
} // namespace Dune

630
#endif // DUNE_DORIE_MODEL_RICHARDS_HH