simulation.hh 5.96 KB
Newer Older
1 2
#ifndef DUNE_DORIE_SIMULATION_HH
#define DUNE_DORIE_SIMULATION_HH
Lukas Riedel's avatar
Lukas Riedel committed
3

4 5 6 7 8 9 10
#include <memory>
#include <string>
#include <iostream>

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

Santiago Ospina's avatar
Santiago Ospina committed
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/gridoperator/gridoperator.hh>
#include <dune/pdelab/gridoperator/onestep.hh>
#include <dune/pdelab/newton/newton.hh>

#include "util.hh"
#include "output.hh"
#include "adaptivity.hh"
#include "../solver/util_interpolator.hh"
#include "../solver/param_factory.hh"
#include "../solver/richards_boundary.hh"
#include "../solver/richards_source.hh"
#include "../solver/richards_initial.hh"
#include "../solver/operator_DG.hh"
#include "../solver/util_controller.hh"

Lukas Riedel's avatar
Lukas Riedel committed
27 28
namespace Dune{
namespace Dorie{
29

30
/// Basic Simulation class providing objects and functions for computing the solution
31
template<typename Traits>
32
class Simulation
33
{
34
protected:
35 36 37

	static constexpr int dim = Traits::dim;
	static constexpr int order = Traits::fem_order;
38
	using RF = typename Traits::RF;
39 40
	using Grid = typename Traits::Grid;
	using GV = typename Traits::GV;
41 42

	/// GFS Helper
43
	using GFSHelper = Dune::Dorie::GridFunctionSpaceHelper<GV,RF,Traits::fem_order,Traits::GridGeometryType>;
44
	/// Problem GFS
45
	using GFS = typename GFSHelper::Type;
46
	/// GFS Constraints Container
47
	using CC = typename GFSHelper::CC;
48

49 50

	/// LSGFS Helper
51
	using LSGFSHelper = Dune::Dorie::LinearSolverGridFunctionSpaceHelper<GV,RF,Traits::GridGeometryType>;
52 53 54 55 56
	/// Linear solver GFS
	using LSGFS = typename LSGFSHelper::Type;
	/// LSGFS Constraints Container
	using LSCC = typename LSGFSHelper::CC;

57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
	// -- DORiE Class Definitions -- //
	/// Parameter Interpolator
	using Interpolator = Dune::Dorie::InterpolatorBase<Traits,Traits::dim>;
	/// Factory class for Interpolator
	using InterpolatorFactory = Dune::Dorie::InterpolatorFactory<Traits,Traits::dim>;
	/// Class handling soil parameters
	using Parameters = Dune::Dorie::ParameterizationBase<Traits,Interpolator>;
	/// Factory class for Parmeter class
	using ParameterizationFactory = Dune::Dorie::ParameterizationFactory<Traits,Interpolator>;
	/// Class defining boundary conditions
	using FlowBoundary = Dune::Dorie::FlowBoundary<Traits>;
	/// Class defining source terms
	using FlowSource = Dune::Dorie::FlowSource<Traits,Parameters>;
	/// Class defining the initial condition
	using FlowInitial = Dune::Dorie::FlowInitial<Traits,Parameters>;
	/// Local spatial operator
73
	using SLOP = Dune::Dorie::Operator::RichardsDGSpatialOperator<Traits,Parameters,FlowBoundary,FlowSource,typename GFSHelper::FEM,false>;
74
	/// Local temporal operator
75
	using TLOP = Dune::Dorie::Operator::RichardsDGTemporalOperator<Traits,Parameters,typename GFSHelper::FEM,false>;
76 77 78 79 80
	/// Time controller
	using CalculationController = Dune::Dorie::CalculationController<RF,FlowBoundary>;

	// -- DUNE Class Defintions -- //
	/// Matrix backend type
Lukas Riedel's avatar
Lukas Riedel committed
81
	using MBE = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
82
	/// Grid operator for spatial LOP
83
	using GO0 = Dune::PDELab::GridOperator<GFS,GFS,SLOP,MBE,RF,RF,RF,CC,CC>;
84
	/// Grid operator for temporal LOP
85
	using GO1 = Dune::PDELab::GridOperator<GFS,GFS,TLOP,MBE,RF,RF,RF,CC,CC>;
86 87 88 89
	/// One step grid operator for merging both GOs
	using IGO = Dune::PDELab::OneStepGridOperator<GO0,GO1>;
	/// Solution vector type
	using U = typename IGO::Traits::Domain;
90
	/// Linear solver types
91
	using LS = Dune::PDELab::ISTLBackend_OVLP_AMG_4_DG<IGO,CC,LSGFS,LSCC,
92
		Dune::PDELab::CG2DGProlongation,Dune::SeqSSOR,Dune::BiCGSTABSolver>;
93
	/// Non-linear solver types
94
	using PDESOLVER = Dune::PDELab::Newton<IGO,LS,U>;
95 96
	/// Time stepping scheme
	using TimeStepScheme = Dune::PDELab::Alexander2Parameter<RF>;
97
	/// Methods computing the time step
98
	using OSM = Dune::PDELab::OneStepMethod<RF,IGO,PDESOLVER,U,U>;
99

100 101
	// -- Utility Class Definitions -- //
	/// VTK Output writer base class
102
	using OutputWriter = Dune::Dorie::OutputWriterBase<Traits>;
103
	/// Factory for creating the output writer
104
	using OutputWriterFactory = Dune::Dorie::OutputWriterFactory<Traits,GFS,Parameters,U>;
105
	/// Grid Adaptivity Handler base class
106
	using AdaptivityHandler = Dune::Dorie::AdaptivityHandlerBase<Traits,GFS,Parameters,FlowBoundary,U>;
107
	/// Factory for creating the AdaptivityHandler
108
	using AdaptivityHandlerFactory = Dune::Dorie::AdaptivityHandlerFactory<Traits,GFS,Parameters,FlowBoundary,U>;
109

110 111
	Dune::MPIHelper& helper;
	std::shared_ptr<Grid> grid;
112
	Dune::ParameterTree& inifile;
113
	GV gv;
114 115 116
	std::unique_ptr<GFS> gfs;
	std::unique_ptr<CC> cc;

117 118 119
	std::unique_ptr<LSGFS> lsgfs;
	std::unique_ptr<LSCC> lscc;

120 121 122 123 124 125 126 127
	std::unique_ptr<Interpolator> interp;
	std::unique_ptr<Parameters> param;
	std::unique_ptr<FlowBoundary> fboundary;
	std::unique_ptr<FlowSource> fsource;
	std::unique_ptr<FlowInitial> finitial;
	std::unique_ptr<CalculationController> controller;

	std::unique_ptr<SLOP> slop;
Lukas Riedel's avatar
Lukas Riedel committed
128
	MBE mbe_slop;
129
	std::unique_ptr<TLOP> tlop;
Lukas Riedel's avatar
Lukas Riedel committed
130
	MBE mbe_tlop;
131 132 133
	std::unique_ptr<GO0> go0;
	std::unique_ptr<GO1> go1;
	std::unique_ptr<IGO> igo;
134 135 136
	std::unique_ptr<LS> ls;
	std::unique_ptr<PDESOLVER> pdesolver;
	std::unique_ptr<OSM> osm;
137 138 139 140
	TimeStepScheme alex2;

	std::unique_ptr<U> uold;
	std::unique_ptr<U> unew;
141

142 143
	std::unique_ptr<OutputWriter> output;
	std::unique_ptr<AdaptivityHandler> adaptivity;
144

145 146
	const int verbose;

147
public:
148

149
	/// Execute the simulation until tEnd is reached.
150
	void run ();
151

152 153 154 155 156 157
	/// Construct basic simulation. Prepare setup for basic simulation.
	/**
	 *  \param _helper Dune MPI instance controller
	 *  \param _grid Shared pointer to the grid
	 *  \param _inifile Dune parameter file parser
	 */
158
	Simulation (Dune::MPIHelper& _helper, std::shared_ptr<Grid> _grid, Dune::ParameterTree& _inifile);
159 160 161 162

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

165 166 167 168 169
	/// Compute a time step
	/** Catch errors in the Newton and the Linear Solver.
	 *  \throw Dune::Exception Fatal error occurs during computation
	 *  \return True if computation succeeded
	 */
170
	bool compute_time_step ();
171

172 173
};

Lukas Riedel's avatar
Lukas Riedel committed
174 175
} // namespace Dorie
} // namespace Dune
176

177
#endif // DUNE_DORIE_SIMULATION_HH