root/src/secOrd_op_initData_impl.cpp

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
  1. r_div_rEnd
  2. eoc_velocity
  3. Explicit_initial_data
  4. Random_initial_data
  5. Random_initial_data
  6. clone
  7. interpolate
  8. evaluate
  9. clone
  10. interpolate
  11. evaluate
  12. clone
  13. interpolate
  14. evaluate
  15. interpolate
  16. evaluate
  17. interpolate
  18. evaluate
  19. interpolate
  20. evaluate
  21. evaluate
  22. hom_value
  23. pertubation
  24. print_configuration
  25. dof_filename

   1 /*! \file secOrd_op_initData_impl.cpp
   2     \brief Implementation of secOrd_op_initData_impl.h
   3 
   4      Revision history
   5      --------------------------------------------------
   6 
   7           Revised by Christian Power June 2016
   8           Revised by Christian Power April 2016
   9           Originally written by Christian Power
  10                (power22c@gmail.com) Februar 2016
  11 
  12      Idea
  13      --------------------------------------------------
  14 
  15      Implementing `Explicit_initial_data` and `Random_initial_data`.
  16 
  17      \author Christian Power 
  18      \date 15. June 2016
  19      \copyright Copyright (c) 2016 Christian Power.  All rights reserved.
  20  */
  21 
  22 #include <iostream>
  23 #include <sstream>
  24 #include <fstream>
  25 #include <cmath>
  26 #include <numeric>
  27 #include <config.h>
  28 #include <dune/fem/io/parameter.hh>
  29 #include <dune/fem/operator/lagrangeinterpolation.hh>
  30 #include "secOrd_op_initData_impl.h"
  31 #include "io_parameter.h"
  32 #include "esfem_error.h"
  33 
  34 
  35 using Esfem::Impl::Explicit_initial_data;
  36 using Esfem::Impl::Random_initial_data;
  37 using Esfem::Impl::Analytic_velocity;
  38 using Esfem::Impl::sphere_1EF;
  39 using Esfem::Impl::sphere_2EF;
  40 using Esfem::Impl::sphere_3EF;
  41 using Esfem::Impl::sphere_eigenFun;
  42 using Esfem::Impl::sphere_mcf_sol;
  43 using Esfem::Impl::sls_iData;
  44 using Esfem::Impl::sls_v_iData;
  45 using Esfem::Impl::sd_iData;
  46 using Dune::Fem::Parameter;
  47 using namespace std;
  48 //! \f$ \R^3 \f$
  49 using Vec_domain = Analytic_velocity::Domain;
  50 //! \f$ \R^3 \f$
  51 using Vec_range = Analytic_velocity::Range;
  52 
  53 // ----------------------------------------------------------------------
  54 // Some static inline analytic expression
  55 
  56 //! Helper for eoc_velocity()
  57 /*! \param t Current time \f$t\f$
  58   \returns \f$r(t)/r_{end} = \frac{r_0}{r_{end} e^{-kt} + r_0 (1-e^{-kt})}\f$
  59 */
  60 static inline double r_div_rEnd(const double t){
  61   const double r0 = 1., r_end = 2., kt = .5 * t;
  62   return r0 / (r_end * exp(-kt) + r0 * (1 - exp(-kt)) );
  63 }
  64 //! Velocity for the solution driven 2016 paper
  65 /*! Eoc means that this velocity was used to generate the eoc tables.
  66   \param t Current time \f$t\f$
  67   \param d Position  \f$x\f$ at time \f$t\f$
  68   \retval r \f$v(x,t) = k \Bigl(1 - \frac{r(t)}{r_{end}}\Bigr) x\f$
  69   \sa r_div_rEnd
  70  */
  71 static inline void eoc_velocity(const double t, const Vec_domain& d, Vec_range& r){
  72   const double k = .5, // should be read from a parameter file
  73     r_d_rEnd = r_div_rEnd(t),
  74     factor =   k * (1 - r_d_rEnd);
  75   r[0] = d[0] * factor;
  76   r[1] = d[1] * factor;
  77   r[2] = d[2] * factor;
  78 }
  79 
  80 // ----------------------------------------------------------------------
  81 // Implementation Explicit_initial_data
  82 
  83 Explicit_initial_data::
  84 Explicit_initial_data(const Esfem::Grid::Grid_and_time& gt,
  85                       const Esfem::Growth type)
  86   :tp {gt.time_provider()}
  87 {
  88   switch(type){
  89   case Growth::promoting:
  90     fun_impl = // [tp_ptr = &tp]
  91       [&tp = tp](const Domain& d, Range& r){
  92       const double x = d[0];
  93       const double y = d[1];
  94       const double t = tp.time();
  95       r = x * y * std::exp(-6. * t);
  96     };
  97     break;
  98   case Growth::inhibiting:
  99     fun_impl = [&tp = tp](const Domain& d, Range& r){
 100       const double y = d[1];
 101       const double z = d[2];
 102       const double t = tp.time();
 103       r = y * z * std::exp(-6. * t);
 104     };
 105     break;
 106   default:
 107     throw InitData_error {Assert::compose(__FILE__, __LINE__, "Bad Growth type")};
 108     break;
 109   };
 110 }
 111 
 112 // ----------------------------------------------------------------------
 113 // Implementation Random_initial_data
 114 
 115 Random_initial_data::
 116 Random_initial_data(const Esfem::Io::Parameter& p,
 117                     const Esfem::Growth type)
 118   :Random_initial_data {hom_value(p, type), pertubation(p, type)}
 119 {
 120   std::cout << print_configuration(p, type) << std::endl;
 121 }
 122 Random_initial_data::
 123 Random_initial_data(const double hom_value,
 124                     const double pertubation)
 125   :random_fun {std::bind(Random_dist {hom_value, hom_value + pertubation},
 126                          Random_engine {})}
 127 {}
 128 
 129 // ----------------------------------------------------------------------
 130 // sphere_1EF
 131 
 132 sphere_1EF::sphere_1EF(const Grid::Grid_and_time& gt)
 133   :tp {gt.time_provider()} {}
 134 sphere_1EF* sphere_1EF::clone(){
 135   return new sphere_1EF {*this};
 136 }
 137 void sphere_1EF::interpolate(Grid::Scal_FEfun& fef) const{
 138   using Fef = Grid::Scal_FEfun::Dune_FEfun;
 139   Dune::LagrangeInterpolation<Fef>::interpolateFunction(*this, fef); 
 140 }
 141 void sphere_1EF::evaluate(const domain& x, range& y)const{
 142   y = x[0]*x[1]*exp(-6*tp.time());
 143 }
 144 
 145 // ----------------------------------------------------------------------
 146 // sphere_eigenFun
 147 
 148 sphere_eigenFun::sphere_eigenFun(const Grid::Grid_and_time& gt)
 149   :tp {gt.time_provider()} {}
 150 sphere_eigenFun* sphere_eigenFun::clone(){
 151   return new sphere_eigenFun {*this};
 152 }
 153 void sphere_eigenFun::interpolate(Grid::Vec_FEfun& vfef) const{
 154   using Vfef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
 155   Dune::LagrangeInterpolation<Vfef>::interpolateFunction(*this, vfef); 
 156 }
 157 void sphere_eigenFun::evaluate(const Domain& x, Range& y) const{
 158   y[0] = x[0]*x[1]; // xy
 159   y[1] = x[1]*x[2]; // yz
 160   y[2] = x[0]*x[2]; // xz
 161   y *= exp(-6.*tp.time());
 162 }
 163 
 164 // ----------------------------------------------------------------------
 165 // sphere_mcf_sol
 166 
 167 sphere_mcf_sol::sphere_mcf_sol(const Grid::Grid_and_time& gt) 
 168   :tp {gt.time_provider()} {}
 169 sphere_mcf_sol* sphere_mcf_sol::clone(){
 170   return new sphere_mcf_sol {*this};
 171 }
 172 void sphere_mcf_sol::interpolate(Grid::Vec_FEfun& rhs) const{
 173   using vfef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
 174   Dune::LagrangeInterpolation<vfef>::interpolateFunction(*this, rhs);
 175 }
 176 void sphere_mcf_sol::evaluate(const Domain& x, Range& y) const{
 177   auto norm = 0.;
 178   for(int i = 0; i < Domain::dimension; ++i) norm += x[i]*x[i];
 179   norm = sqrt(norm);
 180   const auto rt = sqrt(1. - 2 * 2 * tp.time());
 181   y = x;
 182   y *= rt / norm;
 183 }
 184 
 185 // ----------------------------------------------------------------------
 186 // sls_iData
 187 
 188 sls_iData::sls_iData(const Grid::Grid_and_time& gt)
 189   :tp {gt.time_provider()},
 190    rA {Parameter::getValue<double>("logistic_growth.r_start", 1.)},
 191    rE {Parameter::getValue<double>("logistic_growth.r_end", 2.)},
 192    k {Parameter::getValue<double>("logistic_growth.steepness", .5)}
 193 {}  
 194 void sls_iData::interpolate(Grid::Vec_FEfun& rhs) const{
 195   using vfef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
 196   Dune::LagrangeInterpolation<vfef>::interpolateFunction(*this, rhs);
 197 }
 198 void sls_iData::evaluate(const Domain& x, Range& y) const{
 199   const auto 
 200     norm = sqrt(inner_product(&x[0], &x[0]+Domain::dimension, &x[0], 0.)),
 201     ekt = exp(-k*tp.time()),
 202     lgf = rE*rA/(rE * ekt + rA * (1 - ekt));
 203   y = x;
 204   y *= lgf / norm;
 205 }
 206 
 207 sls_v_iData::sls_v_iData(const Grid::Grid_and_time& gt)
 208   :tp {gt.time_provider()},
 209    rA {Parameter::getValue<double>("logistic_growth.r_start", 1.)},
 210    rE {Parameter::getValue<double>("logistic_growth.r_end", 2.)},
 211    k {Parameter::getValue<double>("logistic_growth.steepness", .5)}
 212 {}  
 213 void sls_v_iData::interpolate(Grid::Vec_FEfun& rhs) const{
 214   using vfef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
 215   Dune::LagrangeInterpolation<vfef>::interpolateFunction(*this, rhs);
 216 }
 217 void sls_v_iData::evaluate(const Domain& x, Range& y) const{
 218   const auto 
 219     ekt = exp(-k*tp.time()),
 220     rt = rE*rA/(rE * ekt + rA * (1. - ekt));
 221   y = x;
 222   y *= k * (1. - rt/ rE);
 223 }
 224 
 225 sd_iData::sd_iData(const Grid::Grid_and_time& gt) :tp {gt.time_provider()} {}
 226 void sd_iData::interpolate(Grid::Vec_FEfun& rhs) const{
 227   using vfef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
 228   Dune::LagrangeInterpolation<vfef>::interpolateFunction(*this, rhs);
 229 }
 230 void sd_iData::evaluate(const Domain& x, Range& y) const{
 231   const auto
 232     norm = sqrt(inner_product(&x[0], &x[0]+Domain::dimension, &x[0], 0.)),
 233     fac = exp(tp.time()) / norm;
 234   y = x;
 235   y *= fac;
 236 }
 237 // ----------------------------------------------------------------------
 238 // Analytic_velocity
 239 
 240 Analytic_velocity::Analytic_velocity(const Esfem::Grid::Grid_and_time& gt)
 241   :tp {gt.time_provider()}
 242 {}
 243 
 244 void Analytic_velocity::evaluate(const Domain& d, Range& r) const{
 245   const double t = tp.time();
 246   eoc_velocity(t, d, r);
 247 }
 248 
 249 // ----------------------------------------------------------------------
 250 // helper functions
 251 
 252 double Esfem::Impl::
 253 hom_value(const Esfem::Io::Parameter& p, const Esfem::Growth type){
 254   double rv = 0.;
 255   switch(type){
 256   case Esfem::Growth::promoting:
 257     rv = p.u_hom_value();
 258     break;
 259   case Esfem::Growth::inhibiting:
 260     rv = p.w_hom_value();
 261     break;
 262   default:
 263     throw InitData_error{Assert::compose(__FILE__, __LINE__, "Bad Growth type")};
 264     break;
 265   };
 266   return rv;
 267 }
 268 
 269 double Esfem::Impl::
 270 pertubation(const Esfem::Io::Parameter& p, const Esfem::Growth type){
 271   double rv = 0.;
 272   switch(type){
 273   case Esfem::Growth::promoting:
 274     rv = p.u_pertubation();
 275     break;
 276   case Esfem::Growth::inhibiting:
 277     rv = p.w_pertubation();
 278     break;
 279   default:
 280     throw InitData_error{Assert::compose(__FILE__, __LINE__, "Bad Growth type")};
 281     break;    
 282   };
 283   return rv;
 284 }
 285 
 286 std::string Esfem::Impl::print_configuration(const Esfem::Io::Parameter& p,
 287                                              const Esfem::Growth type){
 288   std::ostringstream oss;
 289   oss << "Using Random_initial_data():\n"
 290     "Random distribution: uniform_real_distribution<> with\n"
 291       << "hom_value: " << hom_value(p, type) << '\n'
 292       << "pertubation: " << pertubation(p, type);
 293   return oss.str();
 294 }
 295 
 296 std::string Esfem::Impl::dof_filename(const Io::Parameter& p, const Growth type){
 297   std::string rv {};
 298   switch(type){
 299   case Growth::promoting:
 300     rv = p.u_init_dof();
 301     break;
 302   case Growth::inhibiting:
 303     rv = p.w_init_dof();
 304     break;
 305   default:
 306     throw InitData_error {Assert::compose(__FILE__, __LINE__, "Bad Growth type")};
 307     break;
 308   };
 309   return rv;
 310 }
 311 
 312 // ----------------------------------------------------------------------
 313 // Implementation Init_data::Data
 314 
 315 Esfem::SecOrd_op::Init_data::Data::Data(const Grid::Grid_and_time& gt,
 316                                         const Growth type)
 317   :eid_ptr {std::make_unique<Explicit_initial_data>(gt,type)}
 318 {}
 319 
 320 Esfem::SecOrd_op::Init_data::Data::Data(const Io::Parameter& p, const Growth type)
 321   :dof_io_filename {Impl::dof_filename(p, type)},
 322    rid_ptr {std::make_unique<Random_initial_data>(p, type)}
 323 {}
 324 
 325 // ----------------------------------------------------------------------
 326 // Exact_velocity::Data
 327 
 328 Esfem::SecOrd_op::Exact_velocity::Data::Data(const Grid::Grid_and_time& gt)
 329   :v_fun {gt}
 330 {}

/* [<][>][^][v][top][bottom][index][help] */