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