The TS-GitLab will have to shut down towards the end of the year — please think about migrating your projects to GitLab.com or GitHub.
(This is still a very early message, meant to keep you informed. There will be more urgent ones in the future.)

knofu.hh 6.49 KB
Newer Older
Lukas Riedel's avatar
Lukas Riedel committed
1 2 3 4 5 6
#ifndef DUNE_DORIE_KNOFU_HH
#define DUNE_DORIE_KNOFU_HH

namespace Dune{
namespace Dorie{

7 8 9 10 11 12 13 14
template<typename Traits, typename DGF, typename GV, typename Parameters>
class MatricHeadAdapter : public Dune::PDELab::GridFunctionBase
	<Dune::PDELab::GridFunctionTraits<GV,typename Traits::RangeField,1,typename Traits::Scalar>,
		MatricHeadAdapter<Traits,DGF,GV,Parameters> >
	{
	private:
		using GFTraits = Dune::PDELab::GridFunctionTraits<GV,typename Traits::RangeField,1,typename Traits::Scalar>;
	public:
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
			/*
			* \brief Constructor
			* \param _gv GridView
			* \param _dgf grid function to interpolate from
			* \param _param Parameter object
			*/
			MatricHeadAdapter (const GV& _gv, const DGF& _dgf, const Parameters& _param) : 
				gv(_gv), dgf(_dgf), param(_param)
			{}

			/**
			* \brief Evaluation of grid function at certain position
			* \param e Element entity
			* \param x Position in local entity coordinates
			* \param y Return value
			* \return Value of grid function
			*/
			void evaluate (const typename GFTraits::ElementType& e,
				const typename GFTraits::DomainType& x,
				typename GFTraits::RangeType& y) const
			{
				typename GFTraits::RangeType u;
				dgf.evaluate(e,x,u);

				// calculate matric head
				const typename GFTraits::DomainType xGlobal = e.geometry().global(x);
				const typename GFTraits::RangeType head = param.matricHead(u,xGlobal);
				y = param.headRefToMiller(head,xGlobal);
			}

			const GV& getGridView () const
			{
				return gv;
			}

		private:
			const GV& gv;
			const Parameters& param;
			const DGF& dgf;
		};
55

Lukas Riedel's avatar
Lukas Riedel committed
56 57 58
template<typename Traits, typename SimTraits>
struct KnoFuInterfaceTraits
{
59 60 61
private:
	using RF = typename Traits::RF;

Lukas Riedel's avatar
Lukas Riedel committed
62
public:
63
	
64 65 66
};

template<typename Traits>
67
class KnoFuInterface : public RichardsSimulation<Traits>
68
{
69
protected:
70 71 72
	using RF = typename Traits::RF;
	using Grid = typename Traits::Grid;

73
	using Base = RichardsSimulation<Traits>;
74

75 76 77 78 79 80 81 82 83 84 85 86
	using Parameters = typename Base::Parameters;
	/// Interface uses GV on constant level
	using KFGV = typename Traits::Grid::LevelGridView;
	/// Interface GFS implements Problem GFS
	using KFGFSHelper = GridFunctionSpaceHelper<KFGV,RF,Traits::fem_order,Traits::GridGeometryType>;
	using KFGFS = typename KFGFSHelper::Type;
	/// Solution vector type
	using KFU = Dune::PDELab::Backend::Vector<KFGFS,RF>;
	/// Discrete Grid Function of KnoFu solution
	using KFUDGF = Dune::PDELab::DiscreteGridFunction<KFGFS,KFU>;
	/// Adapter for translating water content back into matric head
	using HeadDGF = Dune::Dorie::MatricHeadAdapter<Traits,KFUDGF,KFGV,Parameters>; 
Santiago Ospina's avatar
Santiago Ospina committed
87 88 89 90 91

	/// Matric head
	using GFMatricHead = typename Base::GFMatricHead;
	/// Water content
	using GFWaterContent = typename Base::GFWaterContent;
92 93 94 95 96 97 98

	KFGV kfgv; //!< GridView for KnoFuInterface
	KFGFS kfgfs; //!< GFS for KnoFuInterface
	KFU kfu; //!< Solution for KnoFuInterface

public:

99 100
	/// Export raw solution vector type for external usage
	using SolutionContainer = typename KFU::Container;
101

102
	KnoFuInterface (Dune::MPIHelper& _helper, std::shared_ptr<Grid> _grid, Dune::ParameterTree& _inifile)
103
		: RichardsSimulation<Traits>(_helper,_grid,_inifile),
104
			kfgv(_grid->levelGridView(_inifile.get<int>("grid.initialLevel"))),
105
			kfgfs(KFGFSHelper::create(kfgv)),
106
			kfu(kfgfs,0.0)
107 108 109
	{
		static_assert(!Traits::adapt_grid,"Adaptive Grid Refinement does not work with this interface!");
	}
110 111 112

	void solution_to_water_content ()
	{
Santiago Ospina's avatar
Santiago Ospina committed
113 114
		GFMatricHead udgf(this->gfs,this->unew);
		// GFWaterContent waterdgf(this->gv,*this->param,udgf);
115
		// Dune::PDELab::interpolate(waterdgf,kfgfs,kfu);
116 117 118 119 120 121 122 123 124 125
	}

	void water_content_to_solution ()
	{
		KFUDGF kfudgf(kfgfs,kfu);
		HeadDGF headdgf(kfgv,kfudgf,*this->param);
		Dune::PDELab::interpolate(headdgf,*this->gfs,*this->unew);
		*this->uold = *this->unew;
	}

126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
	SolutionContainer cache_solution ()
	{
		//auto x = this->unew->native();
		return Dune::PDELab::Backend::native(*this->unew);
	}

	SolutionContainer cache_water_content ()
	{
		//auto x = this->unew->native();
		return Dune::PDELab::Backend::native(kfu);
	}

	/// Propagate an initial condition over a certain amount of time
	/** Initial condition is expected as water content and will be transformed
	 *  into matric head. After computing the solution, the matric head state
	 *  is transformed back into water content and returned.
	 *  \param initial_condition Initial state as Dune::BlockVector
	 *  \param t_start Starting time
	 *  \param t_end Finish time
	 *  \return Solution (water content) as Dune::BlockVector
	 */
	SolutionContainer propagate (const SolutionContainer& initial_condition,
		const RF t_start, const RF t_end)
	{
		Dune::PDELab::Backend::native(kfu) = initial_condition;
		water_content_to_solution();

		this->controller->set_time(t_start);
		this->controller->set_t_end(t_end);
		this->run();

		solution_to_water_content();
		SolutionContainer res = Dune::PDELab::Backend::native(kfu);
		return res;
	}

162 163 164 165 166 167
	/// Return observation operator for a single position on the grid
	/** Scalar multiplication of the solution vector with this operator
	 *  will yield the observation.
	 *  \param pos Observation position
	 *  \return Observation vector
	 */
168
	std::vector<RF> observation_operator (const typename Traits::Domain pos)
169
	{
170 171
		SolutionContainer kfu_raw = cache_water_content();
		std::vector<RF> res(kfu_raw.size());
172

173 174
		typedef Dune::PDELab::LocalFunctionSpace<KFGFS> LFS;
		typedef Dune::PDELab::LFSIndexCache<LFS> LFSCache;
175
		typedef typename KFU::template ConstLocalView<LFSCache> XView;
176 177 178
		typedef FiniteElementInterfaceSwitch<
					typename Dune::PDELab::LocalFunctionSpace<KFGFS>::Traits::FiniteElementType
			> FESwitch;
179 180 181 182

		LFS lfs(kfgfs);
		LFSCache lfs_cache(lfs);
		XView x_view(kfu);
183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
		std::vector<typename Traits::Scalar> fem_factors(kfgfs.maxLocalSize()); // FEM coefficient container

		// find entity containing the coordinate
		auto it = kfgv.template begin<0>();
		auto pos_local = pos;
		for( ; it != kfgv.template end<0>(); ++it){
			const auto geo_type = it->type();
			const auto& ref = Dune::ReferenceElements<typename Traits::DF,Traits::dim>::general(geo_type);
			pos_local = it->geometry().local(pos);
			if(ref.checkInside(pos_local))
				break;
		}
		if(it == kfgv.template end<0>())
			DUNE_THROW(Dune::Exception,"Measurement position is not inside grid domain!");

		lfs.bind(*it);
		lfs_cache.update();
		x_view.bind(lfs_cache);
		FESwitch::basis(lfs.finiteElement()).evaluateFunction(pos_local,fem_factors);
		for(std::size_t i = 0; i<fem_factors.size(); ++i){
			const auto index = x_view.cache().containerIndex(i)[0];
			res[index] = fem_factors[i][0];
205
		}
206 207

		return res;
208 209
	}

Lukas Riedel's avatar
Lukas Riedel committed
210 211 212 213 214 215
};

}
}

#endif // DUNE_DORIE_KNOFU_HH