richards.cc 10.4 KB
Newer Older
1
#include <dune/dorie/model/richards/richards.hh>
Lukas Riedel's avatar
Lukas Riedel committed
2

3 4 5 6
namespace Dune{
namespace Dorie{

template<typename Traits>
7 8
ModelRichards<Traits>::ModelRichards (
	const Dune::ParameterTree& _inifile,
9 10
	const GridCreator&   _grid_creator,
	Dune::MPIHelper&     _helper
11
):
12
	ModelBase(log_richards,
13 14
				   _inifile.get<std::string>("output.logLevel"),
				   _helper),
15 16
	helper(_helper),
	inifile(_inifile),
17 18
	output_type(_inifile.get<bool>("output.asciiVtk") ?
		Dune::VTK::OutputType::ascii : Dune::VTK::OutputType::base64),
19
	grid(_grid_creator.grid()),
20
	gv(grid->leafGridView()),
21 22
	mbe_slop(estimate_mbe_entries<typename MBE::size_type>(
				Traits::dim,Traits::GridGeometryType)),
23 24
	mbe_tlop(1),
	time_before(0.0),
25
	dt_before(0.0),
26
	enable_fluxrc(_inifile.get<bool>("fluxReconstruction.enable"))
27
{
28

29
	// --- Output Policy ---
30 31 32
  std::string output_policy_str 
      = inifile.get<std::string>("output.policy");

33 34 35
  if (output_policy_str == "endOfRichardsStep") {
    this->set_policy(OutputPolicy::EndOfRichardsStep);
  } else if (output_policy_str == "none") {
36
    this->set_policy(OutputPolicy::None);
37 38 39 40
  } else {
  	this->_log->error("Invalid output policy '{}'", output_policy_str);
    DUNE_THROW(NotImplemented,"Invalid output policy: " << output_policy_str);
  }
41

42
	// --- Grid Function Space ---
43
	this->_log->trace("Setting up GridFunctionSpace");
44
	gfs = std::make_shared<GFS>(GFSHelper::create(gv));
45 46 47 48
	gfs->name("matric_head");
	cc = std::make_unique<CC>();
	Dune::PDELab::constraints(*gfs,*cc,false);

49
	// --- Create the new parameter class
50
	auto element_map = _grid_creator.element_index_map();
51
	fparam = std::make_shared<FlowParameters>(inifile, grid, element_map);
52

53
	// --- Operator Helper Classes ---
54 55
	fboundary = std::make_shared<const FlowBoundary>(inifile);
	fsource = std::make_shared<const FlowSource>(inifile);
56
	finitial = FlowInitialFactory::create(inifile, gv, this->_log);
57

58 59 60 61
	// Read local operator settings
	namespace OP = Dune::Dorie::Operator;
	const auto settings = OP::read_operator_settings(inifile);

62
	// --- Local Operators ---
63 64
	if constexpr (order>0)
	{
65 66
		this->_log->debug("Setting up local grid operators: DG method");

67 68 69
		const auto upwinding = std::get<OP::RichardsUpwinding>(settings);
		const auto method = std::get<OP::RichardsDGMethod>(settings);
		const auto weights = std::get<OP::RichardsDGWeights>(settings);
70

71
		slop = std::make_unique<SLOP>(inifile, fparam, fboundary, fsource,
72
									  method, upwinding, weights);
73
		tlop = std::make_unique<TLOP>(inifile, fparam);
74 75 76
	}
	else {
		this->_log->debug("Setting up local grid operators: FV method");
77 78
		const auto upwinding = std::get<OP::RichardsUpwinding>(settings);
		slop = std::make_unique<SLOP>(fparam, fboundary, upwinding);
79
		tlop = std::make_unique<TLOP>(fparam);
80
	}
81

82 83
	controller = std::make_unique<CalculationController>(
		inifile, *fboundary, this->_log);
84

85
	// --- Solution Vectors and Initial Condition ---
86 87
	u = std::make_shared<U>(*gfs,.0);
	Dune::PDELab::interpolate(*finitial,*gfs,*u);
88

89
	// --- Writer Setup --- //
Santiago Ospina De Los Ríos's avatar
Santiago Ospina De Los Ríos committed
90
	if (output_policy() != OutputPolicy::None)
91
	{
92 93
		const int subsamling_lvl =
			_inifile.get<int>("output.subsamplingLevel", 0);
94
		const auto output_path = inifile.get<std::string>("output.outputPath");
95 96
		this->_log->debug("Creating a VTK writer with subsampling level: {}",
						  subsamling_lvl);
97 98

		Dorie::create_directories(output_path);
99
		const auto subsamling_intervals = Dune::refinementLevels(subsamling_lvl);
100
		auto subsamling_vtk =
101
			std::make_shared<Dune::SubsamplingVTKWriter<GV>>(gv,
102
															 subsamling_intervals);
103
		vtkwriter = std::make_unique<Writer>(subsamling_vtk,
104
							inifile.get<std::string>("output.fileName"),
105
							output_path,
106 107
							"./");
	}
108

109
	// --- Adaptivity Setup --- //
110 111 112 113 114
	if constexpr (AdaptivityHandlerFactory::enabled)
	{
		AdaptivityHandlerFactory adaptivity_fac(inifile,*grid);
		adaptivity = adaptivity_fac.create();
	}
115

116 117 118
  // --- Operator Setup ---
  operator_setup();
  
119
	this->_log->info("Setup complete");
120 121 122
}

template<typename Traits>
123
void ModelRichards<Traits>::operator_setup()
124
{
125
	auto dof_count = gfs->globalSize();
126 127 128 129
	this->_log->debug("Setting up grid operators and solvers");
	if (helper.size() > 1) {
		this->_log->debug("  DOF of this process: {}", dof_count);
	}
130
	dof_count = gv.comm().sum(dof_count);
131
	this->_log->debug("  Total number of DOF: {}", dof_count);
132

133 134 135 136 137 138 139 140 141 142 143 144 145
	Dune::PDELab::constraints(*this->gfs,*this->cc,false);
	// --- Grid Operators ---
	go0 = std::make_unique<GO0>(*gfs,*cc,*gfs,*cc,*slop,mbe_slop);
	go1 = std::make_unique<GO1>(*gfs,*cc,*gfs,*cc,*tlop,mbe_tlop);
	igo = std::make_unique<IGO>(*go0,*go1);

	// --- Solvers ---
	lsgfs = std::make_unique<LSGFS>(LSGFSHelper::create(gv));
	lscc = std::make_unique<LSCC>();
	ls = std::make_unique<LS>(*igo,*cc,*lsgfs,*lscc,1000,0,true,true);

	pdesolver = std::make_unique<PDESOLVER>(*igo,*ls);
	pdesolver->setParameters(inifile.sub("NewtonParameters"));
146
	pdesolver->setVerbosityLevel(0);
147 148 149 150 151 152 153 154

	// --- Time Step Operators ---
	Dune::PDELab::OneStepMethodResult osm_result;
	if(osm){
		osm_result = osm->result(); // cache old result
	}
	osm = std::make_unique<OSM>(alex2,*igo,*pdesolver);
	osm->setResult(osm_result);
155
	osm->setVerbosityLevel(0);
156 157

	gfs->update();
158 159 160
}

template<typename Traits>
161
void ModelRichards<Traits>::step()
162
{
163
	// long time step log for level 'debug'
164 165 166 167 168
	if (this->_log->should_log(spdlog::level::debug)) {
		this->_log->info("Time Step {}:",
						  osm->result().successful.timesteps);
	}

169 170
	bool step_succeed = false;
	while(not step_succeed)
171
	{
172
		// Obtain time variables
173
		const RF time = controller->getTime();
174
		const RF dt = controller->getDT();
175 176 177 178
		bool exception = false;

		try
		{
179
			// long time step log for level 'debug'
180 181 182 183
			if (this->_log->should_log(spdlog::level::debug)) {
				this->_log->debug("  Time {:.2e} + {:.2e} -> {:.2e}",
								  time, dt, time+dt);
			}
184

185 186 187 188
			const auto iter = controller->getIterations();
			_log->trace("  Allowed iterations for Newton solver: {}", iter);
			pdesolver->setMaxIterations(iter);

189
			// create a DOF vector for next step
190 191
			std::shared_ptr<U> unext = std::make_shared<U>(*u);

192
			dwarn.push(false);
193
			osm->apply(time, dt, *u, *unext);
194 195
			dwarn.pop();

196
			// short time step log for level 'info' (after success)
197 198 199 200 201 202 203
			if (not this->_log->should_log(spdlog::level::debug)) {
				this->_log->info("Time Step {}: {:.2e} + {:.2e} -> {:.2e}",
								 osm->result().successful.timesteps-1,
								 time,
								 dt,
								 time+dt);
			}
204

205
			auto result = osm->result();
206
			this->_log->trace("  Matrix assembly: {:.2e}s (total), {:.2e}s (success)",
207 208
							  result.total.assembler_time,
							  result.successful.assembler_time);
209
			this->_log->trace("  Linear solver:   {:.2e}s (total), {:.2e}s (success)",
210 211
							  result.total.linear_solver_time,
							  result.successful.linear_solver_time);
212

Santiago Ospina De Los Ríos's avatar
Santiago Ospina De Los Ríos committed
213
			u = unext;
214 215
		}
		catch (Dune::PDELab::NewtonError &e){
216
			this->_log->debug("  Solver not converged");
217
			this->_log->trace("  Error in Newton solver: {}", e.what());
218 219 220
			exception = true;
		}
		catch (Dune::ISTLError &e){
221
			this->_log->debug("  Solver not converged");
222
			this->_log->trace("  Error in linear solver: {}", e.what());
223 224 225
			exception = true;
		}
		catch (Dune::Exception &e){
226
			this->_log->debug("  Solver not converged");
227
			this->_log->error("  Unexpected error in solver: {}", e.what());
228
			DUNE_THROW(Dune::Exception, "Critical error in Newton solver");
229 230
		}

231
		// saving last times for adaptivity
232
		time_before = time;
233
		dt_before = dt;
234

235
		// controller reacts to outcome of solution
236
		step_succeed = controller->validate(exception);
237
	}
238

239 240
	// invalidate cache for water flux reconstruction
	cache_water_flux_gf_rt.reset();
241

242
	if (this->output_policy() == OutputPolicy::EndOfRichardsStep)
243
		write_data();
244
}
245

246
template<typename Traits>
247
void ModelRichards<Traits>::mark_grid()
248
{
249 250 251
  if (adaptivity_policy() != AdaptivityPolicy::WaterFlux)
    DUNE_THROW(Dune::NotImplemented,"Not known adaptivity policy for transport solver");
	
252 253
	adaptivity->mark_grid(*grid, gv, *gfs, *fparam,
		*fboundary, time_before+dt_before, *u);
254
}
255

256
template<typename Traits>
257
void ModelRichards<Traits>::adapt_grid()
258
{
259 260
	if constexpr (AdaptivityHandlerFactory::enabled)
		adaptivity->adapt_grid(*grid, *gfs, *u);
261 262 263 264
	else {
		this->_log->error("The grid type does not support adaptation");
		DUNE_THROW(NotImplemented, "Grid does not support adaptation");
	}
265 266 267

	// invalidate cache for water flux reconstruction
	cache_water_flux_gf_rt.reset();
268
}
269

270
template<typename Traits>
271
void ModelRichards<Traits>::post_adapt_grid()
272
{
273
	operator_setup();
274
}
275

276
template<typename Traits>
277
void ModelRichards<Traits>::write_data () const
278
{
Santiago Ospina De Los Ríos's avatar
Santiago Ospina De Los Ríos committed
279
	if (vtkwriter)
280
	{
281 282
		// create adapters
		auto head = get_matric_head();
283 284 285 286
		auto wflux = std::make_shared<const GFWaterFlux>(gfs, u, fparam);
		auto cond = std::make_shared<const GFConductivity>(gv, fparam);
		auto wc =  std::make_shared<const GFWaterContent>(head, gv, fparam);
		auto sat = std::make_shared<const GFSaturation>(head, gv, fparam);
287

288
		if (inifile.get<bool>("output.vertexData")) {
289 290 291 292
			vtkwriter->addVertexData(head,"head");
			vtkwriter->addVertexData(cond,"K_0");
			vtkwriter->addVertexData(wc,"theta_w");
			vtkwriter->addVertexData(sat,"Theta");
293 294 295 296

			if constexpr (order > 0) {
				vtkwriter->addVertexData(wflux, "flux");
			}
297
			if constexpr (enable_rt_engine)
298
				if (enable_fluxrc) {
299
					auto wfluxr = get_water_flux_reconstructed();
300
					auto RT_name = "flux_RT" + std::to_string(flux_order);
301
					vtkwriter->addVertexData(wfluxr,RT_name);
302
				}
303 304 305
		}
		// cell data
		else {
306 307 308 309
			vtkwriter->addCellData(head,"head");
			vtkwriter->addCellData(cond,"K_0");
			vtkwriter->addCellData(wc,"theta_w");
			vtkwriter->addCellData(sat,"Theta");
310 311 312 313

			if constexpr (order > 0) {
				vtkwriter->addCellData(wflux, "flux");
			}
314
			if constexpr (enable_rt_engine)
315
				if (enable_fluxrc) {
316
					auto wfluxr = get_water_flux_reconstructed();
317
					auto RT_name = "flux_RT" + std::to_string(flux_order);
318
					vtkwriter->addCellData(wfluxr,RT_name);
319
				}
320
		}
321

322
		try{
323 324 325
			const auto time = controller->getTime();
			this->_log->trace("Writing solution at time {:.2e}", time);
			vtkwriter->write(time, output_type);
326
		}
327 328 329 330
		catch (Dune::Exception& e) {
			this->_log->error("Writing VTK output failed: {}", e.what());
			DUNE_THROW(Dune::IOError, "Cannot write VTK output!");
		}
331
		catch(...){
332 333
			this->_log->error("Writing VTK output failed for unknown reason");
			DUNE_THROW(Dune::IOError, "Cannot write VTK output!");
334 335
		}
		vtkwriter->clear();
336
	} else {
337 338
		this->_log->error("Calling 'write_data' on object without VTK writer");
		DUNE_THROW(Dune::InvalidStateException, "No vtk writer configured!");
339
	}
340 341
}

342
} // namespace Dorie
343
} // namespace Dune