richards.cc 8.39 KB
Newer Older
1
// -*- tab-width: 2; indent-tabs-mode: nil -*-
Dion Haefner's avatar
Dion Haefner committed
2

3 4 5
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
Santiago Ospina's avatar
Santiago Ospina committed
6

7 8 9
#include <string>
#include <iostream>
#include <cstdio>
10
#include <exception>
11 12 13
#include <unistd.h>
#include <sys/stat.h>
#include <sys/types.h>
14 15 16 17 18 19

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

Santiago Ospina's avatar
Santiago Ospina committed
20
#include <dune/grid/yaspgrid.hh>
21
#include <dune/grid/uggrid.hh>
Santiago Ospina's avatar
Santiago Ospina committed
22

23
#include <dune/dorie/common/setup.hh>
24
#include <dune/dorie/common/logging.hh>
25 26
#include <dune/dorie/common/grid_creator.hh>
#include <dune/dorie/model/richards/richards.hh>
Dion Haefner's avatar
Dion Haefner committed
27

28

Dion Haefner's avatar
Dion Haefner committed
29 30 31 32
//===============================================================
// Main program with grid setup
//===============================================================

33
template<typename Traits>
34
using Sim = Dune::Dorie::RichardsSimulation<Traits>;
35 36

template<int dim, int order>
37
using Simplex = Dune::Dorie::RichardsSimulationTraits<Dune::Dorie::BaseTraits<Dune::UGGrid<dim>,
38
	Dune::GeometryType::BasicType::simplex>,order>;
39 40

template<int dim, int order>
41
using Cube = Dune::Dorie::RichardsSimulationTraits<Dune::Dorie::BaseTraits<Dune::YaspGrid<dim>,
42
	Dune::GeometryType::BasicType::cube>,order>;
43 44

template<int dim, int order>
45
using CubeAdaptive = Dune::Dorie::RichardsSimulationTraits<Dune::Dorie::BaseTraits<Dune::UGGrid<dim>,
46
	Dune::GeometryType::BasicType::cube>,order>;
47

48
/// Main Program Function: Read config file, build grid and call Richards Solver
Dion Haefner's avatar
Dion Haefner committed
49 50 51 52 53
int main(int argc, char** argv)
{
  try{
		Dune::Timer timer;

54 55 56
		// Initialize all the stuff!
		const std::string grt = "Starting DORiE";
		auto [inifile, log, helper] = Dune::Dorie::Setup::init(argc, argv, grt);
Dion Haefner's avatar
Dion Haefner committed
57

58 59 60
		// halt process for debugger
		if (inifile.get<bool>("misc.debugMode")) {
			Dune::Dorie::Setup::debug_hook();
61
		}
Dion Haefner's avatar
Dion Haefner committed
62

63
		// setup richards configuration
Santiago Ospina De Los Ríos's avatar
Santiago Ospina De Los Ríos committed
64
		Dune::ParameterTree richards_config = Dune::Dorie::Setup::prep_ini_for_richards(inifile);
65

Dion Haefner's avatar
Dion Haefner committed
66 67 68
		// Read necessary variables
		const std::string gtype = inifile.get<std::string>("grid.gridType");
		const int dim = inifile.get<int>("grid.dimensions");
69
		const std::string adapt_policy_str = inifile.get<std::string>("adaptivity.policy");
70 71
		const int FEorder = richards_config.get<int>("numerics.FEorder");
		const std::string outputPath = richards_config.get<std::string>("output.outputPath");
Dion Haefner's avatar
Dion Haefner committed
72

73 74 75
		Dune::Dorie::AdaptivityPolicy adapt_policy = Dune::Dorie::AdaptivityPolicy::None;

		if (adapt_policy_str == "waterFlux")
76
			adapt_policy = Dune::Dorie::AdaptivityPolicy::WaterFlux;
77
		else if (adapt_policy_str != "none")
78
			DUNE_THROW(Dune::NotImplemented,"The adaptivity policy " <<
79
				adapt_policy_str << " is not implemented!");
80

81

Dion Haefner's avatar
Dion Haefner committed
82
		// Attempt to create output directory
83
		log->info("Creating output directory: {}", outputPath);
Dion Haefner's avatar
Dion Haefner committed
84
		mkdir(outputPath.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
85 86 87
		int result = access(outputPath.c_str(), W_OK);
		if (result != 0)
			DUNE_THROW(Dune::IOError,"Output folder " << outputPath << " not writable");
Dion Haefner's avatar
Dion Haefner committed
88 89 90 91

		if (dim==2)
		{
			if (gtype == "gmsh"){
92
				Dune::Dorie::GridCreator<Dune::UGGrid<2>> grid_creator(inifile, helper);
93 94
				switch(FEorder){
					case 1:{
95
						Sim<Simplex<2,1>> sim(richards_config, grid_creator, helper);
96 97 98
						sim.set_policy(adapt_policy);
						sim.run();
						break;
99
					}
100
					case 2:{
101
						Sim<Simplex<2,2>> sim(richards_config, grid_creator, helper);
102 103 104 105
						sim.set_policy(adapt_policy);
						sim.run();
						break;
					}
106
					case 3:{ // No flux reconstruction available
107
						Sim<Simplex<2,3>> sim(richards_config, grid_creator, helper);
108 109 110
						sim.set_policy(adapt_policy);
						sim.run();
						break;
111
					}
112
					default:
Lukas Riedel's avatar
Lukas Riedel committed
113
						DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
Dion Haefner's avatar
Dion Haefner committed
114 115 116
				}
			}
			else if (gtype == "rectangular"){
117
				if (adapt_policy != Dune::Dorie::AdaptivityPolicy::None) {
118
					Dune::Dorie::GridCreator<Dune::UGGrid<2>> grid_creator(inifile, helper);
Dion Haefner's avatar
Dion Haefner committed
119
					switch(FEorder){
120
						case 1:{
121
							Sim<CubeAdaptive<2,1>> sim(richards_config, grid_creator, helper);
122
							sim.set_policy(adapt_policy);
123
							sim.run();
Dion Haefner's avatar
Dion Haefner committed
124
							break;
125 126
						}
						case 2:{
127
							Sim<CubeAdaptive<2,2>> sim(richards_config, grid_creator, helper);
128
							sim.set_policy(adapt_policy);
129
							sim.run();
Dion Haefner's avatar
Dion Haefner committed
130
							break;
131 132
						}
						case 3:{
133
							Sim<CubeAdaptive<2,3>> sim(richards_config, grid_creator, helper);
134
							sim.set_policy(adapt_policy);
135
							sim.run();
136
							break;
137
						}
Dion Haefner's avatar
Dion Haefner committed
138
						default:
Lukas Riedel's avatar
Lukas Riedel committed
139
							DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
Dion Haefner's avatar
Dion Haefner committed
140 141 142
					}
				}
				else{ // no adaptivity
143
					Dune::Dorie::GridCreator<Dune::YaspGrid<2>> grid_creator(inifile, helper);
Dion Haefner's avatar
Dion Haefner committed
144
					switch(FEorder){
145 146 147 148 149
						case 0:{
							Sim<Cube<2,0>> sim(richards_config, grid_creator, helper);
							sim.run();
							break;
						}
150
						case 1:{
151
							Sim<Cube<2,1>> sim(richards_config, grid_creator, helper);
152
							sim.run();
Dion Haefner's avatar
Dion Haefner committed
153
							break;
154 155
						}
						case 2:{
156
							Sim<Cube<2,2>> sim(richards_config, grid_creator, helper);
157
							sim.run();
Dion Haefner's avatar
Dion Haefner committed
158
							break;
159 160
						}
						case 3:{
161
							Sim<Cube<2,3>> sim(richards_config, grid_creator, helper);
162
							sim.run();
163
							break;
164
						}
Dion Haefner's avatar
Dion Haefner committed
165
						default:
Lukas Riedel's avatar
Lukas Riedel committed
166
							DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
Dion Haefner's avatar
Dion Haefner committed
167 168 169 170 171 172 173 174 175 176
					}
				}
			}
			else
				DUNE_THROW(Dune::NotImplemented,"Grid Type not supported!");
		}

		else if (dim==3)
		{
			if (gtype == "gmsh"){
177
				Dune::Dorie::GridCreator<Dune::UGGrid<3>> grid_creator(inifile, helper);
178 179
				switch(FEorder){
					case 1:{
180
						Sim<Simplex<3,1>> sim(richards_config, grid_creator, helper);
181 182 183
						sim.set_policy(adapt_policy);
						sim.run();
						break;
184
					}
185
					case 2:{
186
						Sim<Simplex<3,2>> sim(richards_config, grid_creator, helper);
187 188 189 190 191
						sim.set_policy(adapt_policy);
						sim.run();
						break;
					}
					case 3:{
192
						Sim<Simplex<3,3>> sim(richards_config, grid_creator, helper);
193 194 195
						sim.set_policy(adapt_policy);
						sim.run();
						break;
196
					}
197
					default:
Lukas Riedel's avatar
Lukas Riedel committed
198
						DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
Dion Haefner's avatar
Dion Haefner committed
199 200 201
				}
			}
			else if (gtype == "rectangular"){
202
				if(adapt_policy != Dune::Dorie::AdaptivityPolicy::None){
203
					Dune::Dorie::GridCreator<Dune::UGGrid<3>> grid_creator(inifile, helper);
Dion Haefner's avatar
Dion Haefner committed
204
					switch(FEorder){
205
						case 1:{
206
							Sim<CubeAdaptive<3,1>> sim(richards_config, grid_creator, helper);
207
							sim.set_policy(adapt_policy);
208
							sim.run();
Dion Haefner's avatar
Dion Haefner committed
209
							break;
210 211
						}
						case 2:{
212
							Sim<CubeAdaptive<3,2>> sim(richards_config, grid_creator, helper);
213
							sim.set_policy(adapt_policy);
214
							sim.run();
Dion Haefner's avatar
Dion Haefner committed
215
							break;
216
						}
217
						case 3:{ // No flux reconstruction available
218
							Sim<CubeAdaptive<3,3>> sim(richards_config, grid_creator, helper);
219
							sim.set_policy(adapt_policy);
220
							sim.run();
221
							break;
222
						}
Dion Haefner's avatar
Dion Haefner committed
223
						default:
Lukas Riedel's avatar
Lukas Riedel committed
224
							DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
Dion Haefner's avatar
Dion Haefner committed
225 226 227
					}
				}
				else{ // no adaptivity
228
					Dune::Dorie::GridCreator<Dune::YaspGrid<3>> grid_creator(inifile, helper);
Dion Haefner's avatar
Dion Haefner committed
229
					switch(FEorder){
230 231 232 233 234
						case 0:{
							Sim<Cube<3,0>> sim(richards_config, grid_creator, helper);
							sim.run();
							break;
						}
235
						case 1:{
236
							Sim<Cube<3,1>> sim(richards_config, grid_creator, helper);
237
							sim.run();
Dion Haefner's avatar
Dion Haefner committed
238
							break;
239 240
						}
						case 2:{
241
							Sim<Cube<3,2>> sim(richards_config, grid_creator, helper);
242
							sim.run();
Dion Haefner's avatar
Dion Haefner committed
243
							break;
244
						}
245
						case 3:{ // No flux reconstruction available
246
							Sim<Cube<3,3>> sim(richards_config, grid_creator, helper);
247
							sim.run();
248
							break;
249
						}
Dion Haefner's avatar
Dion Haefner committed
250
						default:
Lukas Riedel's avatar
Lukas Riedel committed
251
							DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
Dion Haefner's avatar
Dion Haefner committed
252 253 254 255
					}
				}
			}
			else
256
				DUNE_THROW(Dune::NotImplemented,"Grid Type not supported!");
Dion Haefner's avatar
Dion Haefner committed
257 258 259 260 261 262 263
		}

		// grid_dim != 2,3
		else{
		  DUNE_THROW(Dune::NotImplemented,"Number of dimensions (grid.dimensions) not supported!");
		}

264
		log->info("DORiE finished after {:.2e}s :)",
265
					 timer.elapsed());
Dion Haefner's avatar
Dion Haefner committed
266 267

		return 0;
268
	}
269
	catch (Dune::Exception &exc){
270
	  	auto log = Dune::Dorie::create_base_logger();
271
		log->critical("Aborting DORiE after exception: {}", exc.what());
272 273
		return 1;
	}
274
	catch (std::exception& exc) {
275
		auto log = Dune::Dorie::create_base_logger();
276
		log->critical("Aborting DORiE after exception: {}", exc.what());
Dion Haefner's avatar
Dion Haefner committed
277
		return 1;
278 279 280 281
	}
	catch (...){
		auto log = Dune::Dorie::create_base_logger();
		log->error("Unknown error occurred");
282
		log->critical("Aborting DORiE after uncaught exception");
Dion Haefner's avatar
Dion Haefner committed
283
		return 1;
284
	}
Dion Haefner's avatar
Dion Haefner committed
285
}