util_controller.hh 6.27 KB
Newer Older
Dion Haefner's avatar
Dion Haefner committed
1 2 3
#ifndef DUNE_DORIE_CALCULATION_CONTROLLER_HH
#define DUNE_DORIE_CALCULATION_CONTROLLER_HH

4 5 6 7 8 9
#include <iostream>
#include <cmath>

#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>

Dion Haefner's avatar
Dion Haefner committed
10
namespace Dune{
11 12 13 14 15 16 17 18 19 20 21 22 23
namespace Dorie{

/// Class to control time and maximum number of iterations for the Newton solver
template<typename R, typename TSC>
class CalculationController
{
private:

	R time; 		//!< current time
	R dt; 			//!< current time step
	int it; 		//!< current number of allowed iterations for Newton solver
	int verbose; 	//!< verbosity level
	const R dtmax, dtmin; 	//!< time step limits
24
	const R eps; 			//!< error margin
25
	const R dtinc, dtdec; 	//!< time step mutliplicators
26
	R tEnd; 	//!< simulation time limits
27
	const int itmax, itmin; 	//!< newton iteration limits
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42

	const TSC& tsc; //!< Class with Function getNextTimeStamp()


public:

	/// Read relevant time values, check for consistency and adapt time step to BC change, if necessary
	/** \param tsc_ Class with public function getNextTimeStamp(time)
	 */
	CalculationController(const ParameterTree& config, const TSC& tsc_, const Dune::MPIHelper& helper)
		: time(config.get<R>("time.start"))
		, dt(config.get<R>("time.startTimestep"))
		, verbose(config.get<int>("output.verbose"))
		, dtmax(config.get<R>("time.maxTimestep"))
		, dtmin(config.get<R>("time.minTimestep"))
43
		, eps(dtmin/10.0)
44 45 46
		, dtinc(config.get<R>("time.timestepIncreaseFactor"))
		, dtdec(config.get<R>("time.timestepDecreaseFactor"))
		, tEnd(config.get<R>("time.end"))
47
		, itmax(config.get<int>("time.maxIterations"))
48 49 50 51 52 53 54 55 56
		, itmin(config.get<int>("time.minIterations"))
		, tsc(tsc_)
	{
		if(helper.rank()!=0){
			verbose=0;
		}
		userInputCheck();

		// Check whether first time step needs adjustment
57 58
		const R changeTime = tsc.getNextTimeStamp(time);
		if(changeTime>0.0 && time+dt >= changeTime){
59
			dt = changeTime-time-eps;
60 61 62 63 64 65 66 67 68 69 70 71 72 73
			if(verbose>2){
			std::cout 	<< "ADAPTING TIME STEP:  DT : "
						<< std::setprecision(4) << std::scientific << dt
						<< std::endl;
			std::cout 	<< "  Boundary condition change:"
						<< std::setprecision(4) << std::scientific
						<< changeTime << std::endl;
			}
		}

		it=calcIterations(dt);
	}

	/// Return current time
74
	inline R getTime () const { return time; }
75
	/// Return timestep for next calculation
76
	inline R getDT () const { return dt; }
77
	/// Return maximum number of Newton iterations allowed for next calculation
78
	inline int getIterations () const { return it; }
79 80
	/// Return boolean whether new time step shall be computed
	bool doStep () const {
81 82
		if(time<tEnd-eps)
			return true;
83 84
		return false;
	}
Dion Haefner's avatar
Dion Haefner committed
85

86 87 88 89 90 91 92 93 94 95 96
	void set_time (const R _time)
	{
		time = _time;
	}

	/// Set a new finish time for the computation
	void set_t_end (const R _t_end)
	{
		tEnd = _t_end;
	}

97 98 99 100 101 102 103 104 105 106
	/// Validation of solution and time step calculation
	/** If the Newton solver throws an exception (because the solution did not converge
	 *  in the allowed steps or the timestep is too large), exc=true is passed to this
	 *  function and the timestep will be reduced. Returning false causes the Problem
	 *  solver to recalculate the timestep with new parameters provided by the controller.
	 *  The next time step is also adjusted to the next time the BC changes
	 * \param exc Boolean if Newton Class threw an exception
	 * \return Boolean if solution is valid
	 * \throw Dune::Exception Newton solver does not converge at dt>=dtmin and it<=itmax
	 */
107
	bool validate (const bool exc)
108 109 110 111
	{
		bool adapt_to_end = false, adapt_to_bc = false;
		R changeTime = 0.0;

112
		if (!exc) // solution was successfully computed
Dion Haefner's avatar
Dion Haefner committed
113
		{
114
			// time has advanced after calculation
115
			time += dt;
Dion Haefner's avatar
Dion Haefner committed
116

117
			if(time >= tEnd-eps) // Simulation has ended. Do nothing
118 119 120
				return true;

			dt = std::min(dt*dtinc,dtmax);
Dion Haefner's avatar
Dion Haefner committed
121

122
			changeTime = tsc.getNextTimeStamp(time);
123
			if(time >= changeTime-eps && time <= changeTime){
124 125 126 127 128
				// small step to ensure that next time step starts with new BC
				time = changeTime;
				changeTime = tsc.getNextTimeStamp(time);
			}

129
			if(changeTime < 0.0 && time+dt > tEnd){
130 131
				// time step adjustment to tEnd
				dt = tEnd-time;
132
				adapt_to_end = true;
Dion Haefner's avatar
Dion Haefner committed
133
			}
134
			else if(changeTime > 0.0 && time+dt >= changeTime){
135 136
				// time step adjustment to BC Change
				dt = changeTime-time-eps;
137 138 139
				adapt_to_bc = true;
			}
		}
Dion Haefner's avatar
Dion Haefner committed
140

141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
		else 	// Newton solver threw error. Time is not advanced
		{
			if (dt <= dtmin)
				DUNE_THROW(Exception,"Newton Solver does not converge at dtmin and itmax!");
			dt = std::max(dt*dtdec,dtmin);
		}

		// Calculate new iterations based on adapted time step
		it=calcIterations(dt);

		// Output
		if(verbose>2){
			std::cout 	<< "ADAPTING TIME STEP:  DT : "
						<< std::setprecision(4) << std::scientific
						<< dt << std::endl;

			if(adapt_to_end)
				std::cout 	<< "  Simulation end: "
							<< std::setprecision(4) << std::scientific
							<< tEnd << std::endl;
			else if(adapt_to_bc)
				std::cout 	<< "  Boundary condition change: "
							<< std::setprecision(4) << std::scientific
							<< changeTime << std::endl;
			else if(exc)
				std::cout 	<< "  Newton Solver did not converge. Time step is decreased." << std::endl;

			if(dt==dtmax)
				std::cout 	<< "  Maximum time step reached." << std::endl;
			else if(dt==dtmin)
				std::cout 	<< "  Minimum time step reached." << std::endl;

			std::cout 	<< "  Allowed iterations for Newton Solver: " << it << std::endl;
		}

		return !exc;
	}
Dion Haefner's avatar
Dion Haefner committed
178 179


180
private:
Dion Haefner's avatar
Dion Haefner committed
181

182 183 184 185
	/// Calculate the allowed number of Newton iterations by logarithmic interpolation between max and min timestep
	/** \param dt Timestep size
	 *  \return Number of allowed Newton steps for inserted dt
	 */
186
	inline int calcIterations (const R dt)
187
	{
188
		return std::round( (itmin-itmax) * std::log10(dt-dtmin+1.0) / std::log10(dtmax-dtmin+1.0) + itmax  );
189
	}
Dion Haefner's avatar
Dion Haefner committed
190

191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206

	/// Some checks if user input makes sense
	/** \throw Dune::IOError Time and iteration values are not consistent
	 */
	void userInputCheck () const
	{
		if (dtinc<=1.0)
			DUNE_THROW(IOError,"timestepIncreaseFactor must be >1.0!");
		if (dtdec>=1.0 || dtdec<=0.0)
			DUNE_THROW(IOError,"timestepDecreaseFactor exceeds (0.0,1.0)");
		if (itmax<itmin)
			DUNE_THROW(IOError,"maxIterations must be >=minIterations!");
		if (dtmax<dtmin)
			DUNE_THROW(IOError,"maxTimestep must be >=minTimestep!");
		if (dt<dtmin || dt>dtmax)
			DUNE_THROW(IOError,"startTimestep exceeds [minTimestep,maxTimestep]!");
Dion Haefner's avatar
Dion Haefner committed
207
	}
208 209 210 211
};

} // namespace Dorie
} // namespace Dune
Dion Haefner's avatar
Dion Haefner committed
212 213

#endif