adaptivity.hh 8.86 KB
Newer Older
1 2 3 4 5 6
#ifndef DUNE_DORIE_ADAPTIVITY_HH
#define DUNE_DORIE_ADAPTIVITY_HH

namespace Dune{
namespace Dorie{

7
/// Adaptivity base class. Does nothing.
8
template<typename Traits, typename GFS, typename Param, typename U>
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
class AdaptivityHandlerBase
{
private:
	using Grid = typename Traits::Grid;
	using GV = typename Traits::GV;

public:
	AdaptivityHandlerBase () = default;
	virtual ~AdaptivityHandlerBase () = default;

	/// Do nothing.
	/** \return Boolean if operators have to be updated
	 */
	virtual bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, Param& param, U& uold, U& unew) { return false; }
};

25
/// Specialized SimulationBase class providing functions and members for adaptive grid refinement
26 27
template<typename Traits, typename GFS, typename Param, typename U>
class AdaptivityHandler : public AdaptivityHandlerBase<Traits,GFS,Param,U>
28
{
Lukas Riedel's avatar
Lukas Riedel committed
29
private:
30 31 32
	static constexpr int order = Traits::fem_order;
	using RF = typename Traits::RF;
	using Grid = typename Traits::Grid;
33
	using GV = typename Traits::GV;
34

35
	/// Helper class for Adaptivity GFS. implements Problem GFS in order 0
36
	using AGFSHelper = GridFunctionSpaceHelper<GV,RF,0,Traits::GridGeometryType>;
37
	/// Adaptivity GFS
38 39 40 41 42 43 44 45 46 47 48
	using AGFS = typename AGFSHelper::Type;
	/// Error estimator local operator
	using ESTLOP = Dune::Dorie::FluxErrorEstimator<Traits,Param,typename AGFSHelper::FEM>;
	/// Empty constraints container
	using NoTrafo = Dune::PDELab::EmptyTransformation;
	/// Solution vector type
	using U0 = Dune::PDELab::Backend::Vector<AGFS,RF>;
	/// Matrix backend type
	using MBE = Dune::PDELab::istl::BCRSMatrixBackend<RF>;
	/// Grid operator for Error LOP
	using ESTGO = Dune::PDELab::GridOperator<GFS,AGFS,ESTLOP,MBE,RF,RF,RF,NoTrafo,NoTrafo>;
49 50 51

private:

52
	const Dune::ParameterTree& inifile; //!< Parameter file parser
53 54 55 56 57 58
	const int maxLevel; //!< Maximum local grid refinement level
	const int minLevel; //!< Minimum local grid refinement level
	const std::string strategy; //!< Refinement strategy
	const double alpha; //!< Refinement factor
	const double beta; //!< Coarsening factor
	const float adaptivityThreshold; //!< Global error threshold below which no refinement is applied
59
	const int verbose; //!< output verbosity of this object
60 61 62

public:

63
	/// Initialize members from config file parameters.
64 65
	/** Disable grid closure for cubic grid.
	 *  \param _inifile Parameter file parser
66
	 *  \param grid Grid to adapt (reference is not saved)
67
	 */
68
	AdaptivityHandler (Dune::ParameterTree& _inifile, Grid& grid)
69
		: AdaptivityHandlerBase<Traits,GFS,Param,U>(),
70 71 72 73 74 75 76 77
		inifile(_inifile),
		maxLevel(inifile.get<int>("adaptivity.maxLevel")),
		minLevel(inifile.get<int>("adaptivity.minLevel")),
		strategy(inifile.get<std::string>("adaptivity.markingStrategy")),
		alpha(inifile.get<float>("adaptivity.refinementFraction")),
		beta(inifile.get<float>("adaptivity.coarseningFraction")),
		adaptivityThreshold(inifile.get<float>("adaptivity.threshold")),
		verbose(inifile.get<int>("output.verbose"))
78 79 80
	{
		if(strategy!="elementFraction" && strategy!="errorFraction" && strategy!="threshold" && strategy!="targetTolerance")
			DUNE_THROW(Dune::IOError,"Valid values of adaptivity.markingStrategy: elementFraction, errorFraction, threshold, targetTolerance");
81

82
		// QkDG space cannot handle triagles arising from closure of adapted grid
83
		if(Traits::GridGeometryType == Dune::GeometryType::BasicType::cube){
84
			grid.setClosureType(Traits::Grid::ClosureType::NONE);
85
		}
86 87
	}

88 89
	/// Adapt the grid according to the estimated flux error and the strategy applied
	/** 
90 91 92 93 94 95
	 *  \param grid Grid to adapt
	 *  \param gv Leaf grid view of the grid
	 *  \param gfs Solution GridFunctionSpace
	 *  \param param Parameter handler
	 *  \param uold Solution vector
	 *  \param unew Solution vector
96 97 98
	 *  \return Boolean if grid operators have to be set up again (true if grid changed)
	 */
	bool adapt_grid (Grid& grid, GV& gv, GFS& gfs, Param& param, U& uold, U& unew) override
99 100 101 102 103 104
	{
		Dune::Timer timer2;
		float t2;
		double maxeta(0.0);
		int multiRefinementCounter = 0;

105
		if(verbose>1){
106 107 108 109 110 111 112
			std::cout << "GRID ADAPTION [" << strategy << "]:" << std::endl;
		}

		do { // Loop for multiple refinements
			Dune::Timer timer3;
			float t_setup,t_est,t_sqrt,t_strtgy,t_mark,t_adapt;

113
			if(verbose>1 /*&& strategy == "targetTolerance"*/){
114 115 116 117 118 119 120
				std::cout << "  Refinement Step " << multiRefinementCounter << ":  ";
			}

			double eta_alpha(0.0);
			double eta_beta(0.0);

			// Compute error estimate eta
121
			AGFS p0gfs = AGFSHelper::create(gv);
122
			ESTLOP estlop(inifile,param);
123
			MBE mbe(0.0);
124
			ESTGO estgo(gfs,p0gfs,estlop,mbe);
125 126 127 128 129
			U0 eta(p0gfs,0.0);

			t_setup = timer3.elapsed();
			timer3.reset();

130
			estgo.residual(uold,eta);
131 132 133 134 135 136 137 138 139 140

			t_est = timer3.elapsed();
			timer3.reset();

			using Dune::PDELab::Backend::native;
			maxeta = native(eta)[0];
			for (unsigned int i=0; i<eta.flatsize(); i++) {
				native(eta)[i] = sqrt(native(eta)[i]); // eta contains squares
				if (native(eta)[i] > maxeta) maxeta = native(eta)[i];
			}
141
			if (verbose>1) {
142 143 144 145 146 147 148 149 150 151
				std::cout << "Max Local Error: " << maxeta << "   " << std::endl;
				if (maxeta < adaptivityThreshold) {
					std::cout << "  Max local error smaller than threshold " << adaptivityThreshold << ". Skipping grid refinement." << std::endl;
				}
			}

			t_sqrt = timer3.elapsed();
			timer3.reset();

			// Stop adaptivity if target error is reached
152
			if ((maxeta < alpha) && strategy == "targetTolerance") {
153
				if(verbose>1)
154
					std::cout << "  Target error " << alpha << " reached!" << std::endl;
155 156 157 158 159
				break;
			}

			// Apply marking strategy
			if (strategy == "elementFraction")
160
				Dune::PDELab::element_fraction( eta, alpha, beta, eta_alpha, eta_beta, verbose-1 );
161
			else if (strategy == "errorFraction")
162
				Dune::PDELab::error_fraction( eta, alpha, beta, eta_alpha, eta_beta, verbose-1 );
163 164 165 166 167 168 169 170 171 172 173
			else{ //((strategy == "threshold") || (strategy == "targetTolerance")) {
				eta_alpha = alpha;
				eta_beta  = beta; }

			t_strtgy = timer3.elapsed();
			timer3.reset();

			// Skip refinement if threshold is met, but still allow coarsening
			if (maxeta < adaptivityThreshold)
				eta_alpha = maxeta + 1; // Only refine elements with error > maxerror + 1

174
			Dune::PDELab::mark_grid( grid, eta, eta_alpha, eta_beta, minLevel, maxLevel, verbose-1);
175 176 177 178

			t_mark = timer3.elapsed();
			timer3.reset();

179
			Dune::PDELab::adapt_grid( grid, gfs, uold, unew, (order) * 2 );
180

181 182 183 184 185
			t_adapt = timer3.elapsed();
			timer3.reset();

			multiRefinementCounter++;

186
			if(this->verbose>3){
187 188 189
				ios_base_all_saver restore(std::cout);
				std::cout << std::setprecision(4) << std::scientific;
				std::cout << "::: setup time" << std::setw(12) << t_setup << std::endl;
190 191 192 193 194 195 196 197 198
				std::cout << "::: error estimation time" << std::setw(12) << t_est << std::endl;
				std::cout << "::: unsquaring errors time" << std::setw(12) << t_sqrt << std::endl;
				std::cout << "::: strategy application time" << std::setw(12) << t_strtgy << std::endl;
				std::cout << "::: grid marking time" << std::setw(12) << t_mark << std::endl;
				std::cout << "::: grid adaptation time" << std::setw(12) << t_adapt << std::endl;
			}

		} while (strategy == "targetTolerance");

199
		if (verbose>1) {
200
			t2 = timer2.elapsed();
201 202
			ios_base_all_saver restore(std::cout);
			std::cout << std::setprecision(4) << std::scientific;
203
			std::cout << "  Largest local error: " << maxeta << std::endl;
204 205
			std::cout << "::: total adaptation time" << std::setw(12) << t2 << std::endl;
		}
206

207 208 209 210 211 212 213 214 215 216
		return true;
	}

};

/// Factory class for creating an AdaptivityHandler based on Traits
/** \tparam Traits Dune::Dorie::BaseTraits defining the Simulation
 *  The constexpr boolean adapt_grid in Dune::Dorie::BaseTraits will
 *  define which create() function is enabled at compile time.
 */
217
template<typename Traits, typename GFS, typename Param, typename U>
218 219 220 221 222 223
struct AdaptivityHandlerFactory
{
private:
	static constexpr bool enabled = Traits::adapt_grid;
	using Grid = typename Traits::Grid;

224 225
	using AHB = AdaptivityHandlerBase<Traits,GFS,Param,U>;

226 227 228 229 230
	Dune::ParameterTree& inifile;
	Grid& grid;

public:

231 232 233 234 235
	/// Create the factory, taking the parameters for building an AdaptivityHandler
	/** Assert that adaptive grid (UG) is used.
	 *  \param _inifile Parameter file parser
	 *  \param _grid Grid to adapt (reference is not saved)
	 */
236 237 238 239 240 241 242 243 244 245 246
	AdaptivityHandlerFactory (Dune::ParameterTree& _inifile, Grid& _grid) :
		inifile(_inifile),
		grid(_grid)
	{
		static_assert(std::is_same<Grid,Dune::UGGrid<2>>::value ||
									std::is_same<Grid,Dune::UGGrid<3>>::value,
									"Adaptivity is only implemented for UGGrid!");
	}

	/// Create a dummy AdaptivityHandler
	template<bool active = enabled>
247
	std::unique_ptr<AHB> create (typename std::enable_if_t<!active,int> = 0)
248
	{
249
		return std::make_unique<AHB>();
250 251
	}

252 253
	/// Create a true AdaptivityHandler
	template<bool active = enabled>
254
	std::unique_ptr<AHB> create (typename std::enable_if_t<active,int> = 0)
255
	{
256 257
		std::unique_ptr<AHB> p;
		p = std::make_unique<AdaptivityHandler<Traits,GFS,Param,U>>(inifile,grid);
258 259
		return p;
	}
260 261 262 263 264 265 266
};

} // namespace Dorie
} // namespace Dune


#endif // DUNE_DORIE_ADAPTIVITY_HH