src/secOrd_op_initData

/* [<][>][^][v][top][bottom][index][help] */
This source file includes following definitions.
r_div_rEnd
eoc_velocity
Explicit_initial_data
Random_initial_data
Random_initial_data
clone
interpolate
evaluate
clone
interpolate
evaluate
clone
interpolate
evaluate
interpolate
evaluate
interpolate
evaluate
interpolate
evaluate
evaluate
hom_value
pertubation
print_configuration
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] */
root/src/secOrd_op_initData_impl.cpp

DEFINITIONS
