1 #ifndef PROBLEM_TRANSPORT_HH
2 #define PROBLEM_TRANSPORT_HH
3 
4 #include <iostream>
5 #include <sstream>
6 #include <memory>
7 
8 #include <dune/common/fvector.hh>
9 
10 #include "problem.hh"
11 
12 /**
13  * \brief Smooth initial data problem:
14     \f$ \sin( 2\pi x\cdot x ) \f$
15  */
16 template< int dimD >
17 class TransportProblemData1
18 : public ProblemData< dimD, 1 >
19 {
20   typedef ProblemData< dimD, 1 > Base;
21 
22   std::string gridFile_;
23 public:
TransportProblemData1(const int problem)24   TransportProblemData1( const int problem )
25   {
26     if( problem == 3 )
27     {
28       gridFile_ = "/dgf/periodic" + std::to_string(dimDomain) + ".dgf";
29     }
30     else
31     {
32       gridFile_ = "/dgf/unitcube" + std::to_string(dimDomain) + "d.dgf";
33     }
34   }
35 
36   const static int dimDomain = Base::dimDomain;
37   const static int dimRange = Base::dimRange;
38 
39   typedef typename Base::DomainType DomainType;
40   typedef typename Base::RangeType RangeType;
41 
42   //! \copydoc ProblemData::initial
initial(const DomainType & x) const43   RangeType initial ( const DomainType &x ) const
44   {
45     return sin( 2 * M_PI * (x*x) );
46   }
47 
48   //! \copydoc ProblemData::endTime
endTime() const49   double endTime () const
50   {
51     return 0.8;
52   }
53 
gridFile(const std::string & path,const int mpiSize) const54   std::string gridFile ( const std::string &path, const int mpiSize ) const
55   {
56     return path + gridFile_;
57   }
58 
boundaryValue(const DomainType & x,double time) const59   RangeType boundaryValue ( const DomainType &x, double time ) const
60   {
61     return initial( x );
62   }
63 
bndType(const DomainType & normal,const DomainType & x,const double time) const64   int bndType( const DomainType &normal, const DomainType &x, const double time) const
65   {
66     return 1;
67   }
68 
saveInterval() const69   double saveInterval () const
70   {
71     return 0.05;
72   }
73 
74   //! \copydoc ProblemData::refineTol
refineTol() const75   double refineTol () const
76   {
77     return 0.1;
78   }
79   //! \copydoc ProblemData::adaptationIndicator
adaptationIndicator(const DomainType & x,double time,const RangeType & uLeft,const RangeType & uRight) const80   double adaptationIndicator ( const DomainType& x, double time,
81                                const RangeType &uLeft, const RangeType &uRight ) const
82   {
83     return std::abs( uLeft[ 0 ] - uRight[ 0 ] );
84   }
85 };
86 
87 /**
88  * \brief Discontinuous initial data problem: characteristic function for
89     \f$ \{ x\colon |x| < \frac{1}{2} \} \f$
90  */
91 template< int dimD >
92 class TransportProblemData2
93 : public ProblemData< dimD, 1 >
94 {
95   typedef ProblemData< dimD, 1  > Base;
96 
97   std::string gridFile_;
98 public:
TransportProblemData2(const int problem)99   TransportProblemData2( const int problem )
100   {
101     if( problem == 4 )
102     {
103       gridFile_ = "/dgf/periodic" + std::to_string(dimDomain) + ".dgf";
104     }
105     else
106     {
107       gridFile_ = "/dgf/unitcube" + std::to_string(dimDomain) + "d.dgf";
108     }
109   }
110 
111   const static int dimDomain = Base::dimDomain;
112   const static int dimRange = Base::dimRange;
113 
114   typedef typename Base::DomainType DomainType;
115   typedef typename Base::RangeType RangeType;
116 
117   //! \copydoc ProblemData::initial
initial(const DomainType & x) const118   RangeType initial ( const DomainType &x ) const
119   {
120     DomainType r(0.5);
121     return ((x-r).two_norm() < 0.25 ? RangeType( 1 ) : RangeType( 0 ) );
122     //return ((x-r).two_norm() < 0.5 ? RangeType( 1 ) : RangeType( 0 ) );
123   }
124 
125   //! \copydoc ProblemData::endTime
endTime() const126   double endTime () const
127   {
128     return 0.8;
129   }
130 
131   //! \copydoc ProblemData::gridFile
gridFile(const std::string & path,const int mpiSize) const132   std::string gridFile ( const std::string &path, const int mpiSize ) const
133   {
134     return path + gridFile_;
135   }
136 
boundaryValue(const DomainType & x,double time) const137   RangeType boundaryValue ( const DomainType &x, double time ) const
138   {
139     return initial( x );
140   }
141 
bndType(const DomainType & normal,const DomainType & x,const double time) const142   int bndType( const DomainType &normal, const DomainType &x, const double time) const
143   {
144     return 1;
145   }
146 
saveInterval() const147   double saveInterval () const
148   {
149     return 0.05;
150   }
151   //! \copydoc ProblemData::refineTol
refineTol() const152   double refineTol () const
153   {
154     return 0.1;
155   }
adaptationIndicator(const DomainType & x,double time,const RangeType & uLeft,const RangeType & uRight) const156   double adaptationIndicator ( const DomainType& x, double time,
157                                const RangeType &uLeft, const RangeType &uRight ) const
158   {
159     return std::abs( uLeft[ 0 ] - uRight[ 0 ] );
160   }
161 };
162 
163 // TransportModel
164 // ----------------
165 
166 /** \class TransportModel
167  *  \brief description of a transport problem
168  *
169  *  This class describes the following transport problem:
170  *  \f{eqnarray*}
171  *  \partial_t c + \nabla \cdot (v c)
172  *    &=& 0 \quad\mbox{in $\Omega \times ]0,T[$}\\
173  *  c &=& g \quad\mbox{on $\Gamma_{\mathrm{in}}$}\\
174  *  c &=& c_0 \quad\mbox{on $\Omega \times \lbrace 0 \rbrace$}
175  *  \f}
176  */
177 template< int dimD >
178 struct TransportModel
179 {
180   typedef ProblemData< dimD,1 > Problem;
181 
182   typedef typename Problem::DomainType DomainType;
183   typedef typename Problem::RangeType RangeType;
184 
185   static const int dimDomain = Problem::dimDomain;
186   static const int dimRange = Problem::dimRange;
187   static const bool hasFlux = true;
188 
189   /** \brief constructor
190    *  \param problem switch between different data settings
191    */
TransportModelTransportModel192   TransportModel ( int problem )
193   {
194     switch( problem )
195     {
196     case 1:
197     case 3:
198       problem_.reset( new TransportProblemData1< dimDomain >( problem ) );
199       break;
200     case 2:
201     case 4:
202       problem_.reset( new TransportProblemData2< dimDomain >( problem ) );
203       break;
204     default:
205       std::cerr << "Problem " << problem << " does not exists." << std::endl;
206       std::abort();
207     }
208 
209     if( problem < 3 )
210     {
211       // set transport velocity
212       velocity_ = DomainType( 1.25 );
213     }
214     else
215     {
216       velocity_ = 0.0;
217       velocity_[ 0 ] = 1.25;
218     }
219   }
220 
221   /** \brief obtain problem */
problemTransportModel222   const Problem &problem () const
223   {
224     return *problem_;
225   }
226 
fixedDtTransportModel227   double fixedDt () const
228   {
229     return -1;
230   }
231 
232   /** \brief obtain the (constant) velocity for the transport problem */
velocityTransportModel233   const DomainType &velocity () const
234   {
235     return velocity_;
236   }
237 
238   /** \brief evaluate the numerical flux on an intersection
239    *
240    *  \param[in]   normal   scaled normal of the intersection
241    *  \param[in]   time     current time
242    *  \param[in]   xGlobal  evaluation point in global coordinates
243    *  \param[in]   uLeft    value of the solution in the inside entity
244    *  \param[in]   uRight   value of the solution in the outside entity
245    *  \param[out]  flux     numercial flux
246    *
247    *  \returns the maximum wave speed
248    */
numericalFluxTransportModel249   double numericalFlux ( const DomainType &normal,
250                          const double time,
251                          const DomainType &xGlobal,
252                          const RangeType &uLeft, const RangeType &uRight,
253                          RangeType &flux ) const
254   {
255     const double upwind = normal * velocity();
256     flux = upwind * (upwind > 0 ? uLeft : uRight);
257     return std::abs( upwind );
258   }
259 
260   /** \brief evaluate the numerical flux on a boundary
261    *
262    *  \param[in]   normal   scaled normal of the boundary
263    *  \param[in]   time     current time
264    *  \param[in]   xGlobal  evaluation point in global coordinates
265    *  \param[in]   uLeft    value of the solution in the inside entity
266    *  \param[out]  flux     numercial flux
267    *
268    *  \returns the maximum wave speed
269    */
boundaryFluxTransportModel270   double boundaryFlux ( const DomainType &normal,
271                         const double time,
272                         const DomainType &xGlobal,
273                         const RangeType& uLeft,
274                         RangeType &flux ) const
275   {
276     // exact solution is u0(x-ta)
277     DomainType x0( xGlobal );
278     x0.axpy( -time, velocity_ );
279     RangeType uRight = problem().boundaryValue( x0, time );
280     return numericalFlux( normal, time, xGlobal, uLeft, uRight, flux );
281   }
282 
283   /** \brief compute adaptation indicator at intersection
284    *
285    *  \param[in]   normal   scaled normal of the intersection
286    *  \param[in]   time     current time
287    *  \param[in]   xGlobal  evaluation point in global coordinates
288    *  \param[in]   uLeft    value of the solution in the inside entity
289    *  \param[in]   uRight   value of the solution in the outside entity
290    *
291    *  \return value of indicator
292    */
indicatorTransportModel293   double indicator ( const DomainType &normal,
294                      const double time,
295                      const DomainType &xGlobal,
296                      const RangeType &uLeft, const RangeType &uRight) const
297   {
298     return problem().adaptationIndicator( xGlobal, time, uLeft, uRight );
299   }
300 
301   /** \brief compute adaptation indicator at boundary
302    *
303    *  \param[in]   normal   scaled normal of the intersection
304    *  \param[in]   time     current time
305    *  \param[in]   xGlobal  evaluation point in global coordinates
306    *  \param[in]   uLeft    value of the solution in the inside entity
307    *
308    *  \return value of indicator
309    */
boundaryIndicatorTransportModel310   double boundaryIndicator ( const DomainType &normal,
311                              const double time,
312                              const DomainType &xGlobal,
313                              const RangeType& uLeft) const
314   {
315     DomainType x0( xGlobal );
316     x0.axpy( -time, velocity_ );
317     return indicator( normal,time,xGlobal, uLeft, problem().boundaryValue(x0,time) );
318   }
319 
320 protected:
TransportModelTransportModel321   TransportModel ( ) : problem_(), velocity_(1) {}
322 private:
323   std::unique_ptr< Problem > problem_;
324   DomainType velocity_;
325 }; // end class TransportModel
326 
327 #endif // PROBLEM_TRANSPORT_HH
328