root/src/brusselator_algo.cpp

DEFINITIONS

1. brusselator_algo
2. Brusselator_scheme
3. standard_esfem
4. eoc_logisticSphere
5. eoc_mcf
6. eoc_sls
7. sd
8. eoc_sdp
9. prePattern_loop
10. intermediate_action
11. pattern_loop
12. final_action
13. update_surface
14. update_scalar_solution
15. error_on_intSurface
16. pre_loop_action
17. rhs_and_solve_SPDE
18. test

   1 /*! \file brusselator_algo.cpp
   2     \brief Implementation for `brusselator_algo.h`
   3 
   4      Revision history
   5      --------------------------------------------------
   6 
   7           Revised by Christian Power July 2016
   8           Revised by Christian Power June 2016
   9           Revised by Christian Power May 2016
  10           Revised by Christian Power April 2016
  11           Revised by Christian Power March 2016
  12           Originally written by Christian Power
  13                (power22c@gmail.com) February 2016
  14 
  15      Idea
  16      --------------------------------------------------
  17      
  18      Implementation of Esfem::brusselator_algo() and
  19      the class `Esfem::Brusselator_scheme`
  20 
  21      \author Christian Power 
  22      \date 13. July 2016
  23      \copyright Copyright (c) 2016 Christian Power.  All rights reserved.
  24 */
  25 
  26 #include <config.h>
  27 #include "brusselator_algo.h"
  28 #include "brusselator_algo_impl.h"
  29 #include "secOrd_op_solutionDriven.h"
  30 #include "esfem_error.h"
  31 
  32 #ifndef PFILE
  33 #error Give a complete path for parameter file in macro variable PFILE.
  34 #endif
  35 
  36 using Bruss_error = Esfem::BrusselatorScheme_error;
  37 /*!< \brief Shorter name for convenience */
  38 using Esfem::Brusselator_scheme;
  39 using Scal_FEfun_set = Esfem::Grid::Scal_FEfun_set;
  40 /*!< \brief Four functions of type \f$ f\colon \R^3 \to \R \f$ */
  41 using Vec_FEfun_set = Esfem::Grid::Vec_FEfun_set;
  42 /*!< \brief Four functions of type \f$ f\colon \R^3 \to \R^3 \f$ */
  43 using Vector_solver = Esfem::SecOrd_op::Solution_driven;
  44 /*!< \brief Solver for the `X` respectively `surface` variable */
  45 
  46 void Esfem::brusselator_algo(int argc, char** argv){
  47   Dune::Fem::MPIManager::initialize(argc, argv);
  48   constexpr auto parameter_file = PFILE;
  49   Brusselator_scheme fem {argc, argv, parameter_file};
  50   // fem.test();
  51   // fem.prePattern_loop();
  52   // fem.intermediate_action(); 
  53   // fem.pattern_loop();
  54   // fem.final_action();
  55   // fem.standard_esfem(); // code works
  56   // fem.eoc_logisticSphere(); // does not work
  57   // fem.eoc_mcf(); // code works
  58   // fem.eoc_sls(); // code works
  59   // fem.sd(); // code works
  60   fem.eoc_sdp(); // code works
  61 }
  62 
  63 // ----------------------------------------------------------------------
  64 // Brusselator_scheme implementation
  65 
  66 // ------------------------------------------------------------
  67 // Brusselator_scheme public 
  68 
  69 /*! \todo Test if folder for tmp file exists. */
  70 Brusselator_scheme::
  71 Brusselator_scheme(int argc, char** argv,
  72                    const std::string& parameter_fname)
  73 try :data {argc, argv, parameter_fname},
  74   io {data},
  75   fix_grid {data},
  76   norm {fix_grid},  
  77   fef {fix_grid}, 
  78   exact {fix_grid, data} // constructor for random initial data
  79 // exact {fix_grid} // constructor for analytic initial data
  80 {
  81   pre_loop_action(); // initialize member fef
  82 }
  83 catch(...){
  84   throw_with_nested(Bruss_error {"Constructor."});    
  85  }
  86 
  87 void Brusselator_scheme::standard_esfem(){
  88   Rhs load_vector {fix_grid};
  89   Scalar_solver solver {data, fix_grid, fef.u, fef.w};
  90 
  91   error_on_intSurface(); // error on t_0
  92   for(long it = 0; it < time_steps(); ++it){
  93     // rhs = M^n u^n
  94     solver.u.mass_matrix(fef.u.fun, fef.u.rhs_les);
  95     solver.w.mass_matrix(fef.w.fun, fef.w.rhs_les);
  96     next_timeStep(); // next surface
  97 
  98     // rhs += M^{n+1} 1
  99     solver.u.add_massMatrixConstOne_to(fef.u.rhs_les);
 100     solver.w.add_massMatrixConstOne_to(fef.w.rhs_les);
 101     // rhs += load_vector
 102     load_vector.u.assemble_and_addScaled_to(fef.u.rhs_les);
 103     load_vector.w.assemble_and_addScaled_to(fef.w.rhs_les);
 104 
 105     solver.u.solve(fef.u.rhs_les, fef.u.fun);
 106     fef.u.app = fef.u.fun;    
 107     solver.w.solve(fef.w.rhs_les, fef.w.fun);
 108     fef.w.app = fef.w.fun;
 109 
 110     exact.u.interpolate(fef.u.exact);
 111     exact.w.interpolate(fef.w.exact);
 112     error_on_intSurface(); // error on t_n+1
 113   }
 114 }
 115 
 116 void Brusselator_scheme::eoc_logisticSphere(){  
 117   io.identity.interpolate(fef.surface.fun);
 118 
 119   for(long it = 0; it < pattern_timeSteps(); ++it){
 120     Grid::Grid_and_time grid {data,
 121         Grid::compose_dgfName(fef.surface.fun.name(), fef.tmpFile_path), 
 122         fix_grid.time_provider().time()};
 123     Grid::Scal_FEfun_set u {fef.u, grid};
 124     Grid::Vec_FEfun_set X {fef.surface, grid};
 125     SecOrd_op::Init_data u_init {grid, Growth::promoting}; // perhaps not efficient
 126     SecOrd_op::Vec_rhs X_loadVector {grid};
 127     SecOrd_op::Solution_driven X_solver {data, grid, u.app};
 128     
 129     // calcul1ate X.fun
 130     u_init.interpolate(u.app);
 131     X.app = fef.surface.fun;
 132     X_solver.brusselator_rhs(X.app, X.rhs_les); 
 133     // X_loadVector.assemble_and_addScaled_to(X.rhs_les);
 134 
 135     X_solver.solve(X.rhs_les, X.fun);
 136 
 137     next_timeStep();
 138     
 139     // save surface
 140     fef.surface.fun = X.fun; // swap would be more efficient
 141     fef.surface.write(io.dgf_handler, fef.tmpFile_path);    
 142 
 143     // calculate error 
 144     io.identity.interpolate(fef.surface.exact);    
 145     io.surface << fix_grid.time_provider().deltaT() << ' '
 146                << norm.l2_err(fef.surface.fun, fef.surface.exact) << ' '
 147                << norm.h1_err(fef.surface.fun, fef.surface.exact) 
 148                << std::endl;
 149     // I've changed order to get a stationary surface
 150     // io.identity.interpolate(fef.surface.fun);
 151     // fef.surface.write(io.dgf_handler, fef.tmpFile_path);
 152 
 153 
 154     
 155     // update_surface(); // calculate exact surface
 156     // update_scalar_solution(); // on the exact surface 
 157     // rhs_and_solve_SPDE(); // also scalar rhs
 158     // update_velocity();    // exact and approximation 
 159     // error_on_intSurface(); // Error on surface(t_n)
 160     // next_timeStep(); // next surface
 161     // Pattern_helper helper {*this};
 162     // helper.finalize_scalarPDE_rhs();
 163     // helper.solve_scalarPDE();
 164     // helper.errors_on_numSurface();
 165     // helper.plot_paraview();
 166   }
 167 }
 168 /*! \pre delta == 0, epsilon == 1, alpha == 0. */
 169 void Brusselator_scheme::eoc_mcf(){
 170   std::unique_ptr<SecOrd_op::vIdata> ex_ptr 
 171   {SecOrd_op::vIdata::new_sms(fix_grid)};
 172   SecOrd_op::Solution_driven X_solver {data, fix_grid, fef.u.app};
 173   ex_ptr->interpolate(fef.surface.fun);
 174   fef.surface.exact = fef.surface.fun;
 175   for(long it = 0; it < pattern_timeSteps(); ++it){
 176     X_solver.brusselator_rhs(fef.surface.fun, fef.surface.rhs_les);
 177     // les_rhs = M_3 u^n
 178     X_solver.solve(fef.surface.rhs_les, fef.surface.fun);
 179     // (M_3 + tau A_3) u^n+1 = les_rhs
 180     next_timeStep();
 181     fix_grid.new_nodes(fef.surface.exact);
 182     exact.X_ptr->interpolate(fef.surface.exact);
 183     io.surface << fix_grid.time_provider().deltaT() << ' '
 184                << norm.l2_err(fef.surface.fun, fef.surface.exact) << ' '
 185                << norm.h1_err(fef.surface.fun, fef.surface.exact) 
 186                << std::endl;
 187     fix_grid.new_nodes(fef.surface.fun);
 188   }
 189 }
 190 
 191 void Brusselator_scheme::eoc_sls(){
 192   using namespace SecOrd_op;
 193   using std::unique_ptr;
 194   unique_ptr<vRhs> vRhs_ptr {vRhs::new_sls(fix_grid)};
 195   unique_ptr<vIdata> ex_ptr {vIdata::new_sls(fix_grid)};
 196   unique_ptr<sIdata> u_ptr {sIdata::new_1ssef(fix_grid)};
 197   Solution_driven X_solver {data, fix_grid, fef.u.app};
 198   ex_ptr->interpolate(fef.surface.fun);
 199   u_ptr->interpolate(fef.u.app);
 200   fef.surface.exact = fef.surface.fun;
 201   for(long it = 0; it < pattern_timeSteps(); ++it){
 202     X_solver.brusselator_rhs(fef.surface.fun, fef.surface.rhs_les);
 203     // les_rhs = (M_3 + alpha A_3) X^n
 204     vRhs_ptr->addScaled_to(fef.surface.rhs_les);
 205     // les_rhs += tau load_vector^n
 206     X_solver.solve(fef.surface.rhs_les, fef.surface.fun);
 207     // (M_3 + (alpha + tau * epsilon) A_3) X^n+1 = les_rhs
 208     next_timeStep();
 209     fix_grid.new_nodes(fef.surface.exact);
 210     ex_ptr->interpolate(fef.surface.exact);
 211     u_ptr->interpolate(fef.u.app);
 212     io.surface << fix_grid.time_provider().deltaT() << ' '
 213                << norm.l2_err(fef.surface.fun, fef.surface.exact) << ' '
 214                << norm.h1_err(fef.surface.fun, fef.surface.exact) 
 215                << std::endl;
 216     fix_grid.new_nodes(fef.surface.fun);
 217   }
 218 }
 219 
 220 void Brusselator_scheme::sd(){
 221   using namespace SecOrd_op;
 222   std::unique_ptr<vRhs> vRhs_ptr {vRhs::new_sd(fix_grid)};
 223   std::unique_ptr<vIdata> ex_ptr {vIdata::new_sd(fix_grid)};
 224   Solution_driven X_solver {data, fix_grid, fef.u.app};
 225   ex_ptr->interpolate(fef.surface.fun);
 226   fef.surface.exact = fef.surface.fun;
 227   for(long it = 0; it < pattern_timeSteps(); ++it){
 228     X_solver.brusselator_rhs(fef.surface.fun, fef.surface.rhs_les);
 229     vRhs_ptr->addScaled_to(fef.surface.rhs_les);
 230     X_solver.solve(fef.surface.rhs_les, fef.surface.fun);
 231     next_timeStep();
 232     fix_grid.new_nodes(fef.surface.exact);
 233     ex_ptr->interpolate(fef.surface.exact); 
 234     io.surface << fix_grid.time_provider().deltaT() << ' '
 235                << norm.l2_err(fef.surface.fun, fef.surface.exact) << ' '
 236                << norm.h1_err(fef.surface.fun, fef.surface.exact) 
 237                << std::endl;
 238     fix_grid.new_nodes(fef.surface.fun);
 239   }
 240 }
 241 
 242 void Brusselator_scheme::eoc_sdp(){
 243   using namespace SecOrd_op;
 244   using std::unique_ptr;
 245   unique_ptr<vRhs> g_load {vRhs::new_sls(fix_grid)};
 246   unique_ptr<sRhs> f_load {sRhs::new_sdp_u(fix_grid)};
 247   unique_ptr<vIdata> X_ex {vIdata::new_sls(fix_grid)};
 248   unique_ptr<vIdata> V_ex {vIdata::new_v_sls(fix_grid)};
 249   unique_ptr<sIdata> u_ex {sIdata::new_1ssef(fix_grid)};
 250   Scalar_solver solver {data, fix_grid, fef.u, fef.w}; // I just need `solver.u`
 251   Solution_driven X_solver {data, fix_grid, fef.u.app};
 252   X_ex->interpolate(fef.surface.fun);
 253   u_ex->interpolate(fef.u.app);
 254   fef.surface.exact = fef.surface.fun;
 255   fef.u.fun = fef.u.app;
 256   const long end_steps = pattern_timeSteps() + (data.last_step() > data.eps() ? 1 : 0 );
 257   for(long it = 0; it < end_steps; ++it){
 258     fef.velocity.rhs_les = fef.surface.fun; // old surface for velocity
 259     X_solver.brusselator_rhs(fef.surface.fun, fef.surface.rhs_les);
 260     g_load->addScaled_to(fef.surface.rhs_les);
 261     X_solver.solve(fef.surface.rhs_les, fef.surface.fun);
 262     fef.velocity.app = fef.surface.fun; // new surface for velocity 
 263     solver.u.mass_matrix(fef.u.fun, fef.u.rhs_les);
 264     calculate_velocity(fef.velocity.app.cbegin(), fef.velocity.app.cend(),
 265                        fef.velocity.rhs_les.cbegin(), fef.velocity.fun.begin());
 266 
 267     // next surface
 268     next_timeStep();
 269     // if(it < end_steps - 2) next_timeStep();
 270     // else fix_grid.next_timeStep(data.last_step());
 271     
 272     fix_grid.new_nodes(fef.surface.fun);    
 273     f_load->addScaled_to(fef.u.rhs_les);
 274     solver.u.solve(fef.u.rhs_les, fef.u.fun);
 275     fef.u.app = fef.u.fun;
 276 
 277     X_ex->interpolate(fef.surface.exact);
 278     fix_grid.new_nodes(fef.surface.exact); 
 279     u_ex->interpolate(fef.u.exact);
 280     V_ex->interpolate(fef.velocity.exact);
 281     print(io.u, fef.u);
 282     print(io.surface, fef.surface);
 283     print(io.velocity, fef.velocity);
 284     fix_grid.new_nodes(fef.surface.fun);
 285   }
 286 }
 287 
 288 // --------------------------------------------------
 289 // Brusselator_scheme loop action
 290 
 291 void Brusselator_scheme::prePattern_loop(){
 292   PrePattern_helper helper {*this};
 293   
 294 
 295   for(long it = 0; it < prePattern_timeSteps(); ++it, next_timeStep()){
 296     helper.rhs();
 297     helper.solve_pde(); // Does this make any sense?
 298     helper.plot_paraview();
 299   }
 300 }
 301 void Brusselator_scheme::intermediate_action(){
 302   const auto uw_path = fef.tmpFile_path + "intermediate_";
 303   switch(prePattern_timeSteps()){
 304   case 0: // heat.starttime == heat.pattern.endtime
 305     fef.u.read(io.dgf_handler, uw_path);
 306     fef.w.read(io.dgf_handler, uw_path);
 307     break;
 308   default:
 309     fef.u.write(io.dgf_handler, uw_path);
 310     fef.w.write(io.dgf_handler, uw_path);
 311     fef.surface.write(io.dgf_handler, fef.tmpFile_path);
 312   };
 313 }
 314 void Brusselator_scheme::pattern_loop(){
 315   using namespace SecOrd_op;
 316   Scalar_solver solver {data, fix_grid, fef.u, fef.w}; 
 317   Solution_driven X_solver {data, fix_grid, fef.u.app};
 318   Esfem::Io::Paraview paraview {data, fix_grid, fef.u.fun, fef.w.fun};
 319   for(long it = 0; it < pattern_timeSteps(); ++it){
 320     X_solver.brusselator_rhs(fef.surface.fun, fef.surface.rhs_les);
 321     X_solver.solve(fef.surface.rhs_les, fef.surface.fun);
 322     solver.u.mass_matrix(fef.u.fun, fef.u.rhs_les);
 323     solver.w.mass_matrix(fef.w.fun, fef.w.rhs_les);
 324 
 325     next_timeStep(); // next surface
 326     fix_grid.new_nodes(fef.surface.fun);    
 327     solver.u.add_massMatrixConstOne_to(fef.u.rhs_les);
 328     solver.w.add_massMatrixConstOne_to(fef.w.rhs_les);
 329     solver.u.solve(fef.u.rhs_les, fef.u.fun);
 330     fef.u.app = fef.u.fun;
 331     solver.w.solve(fef.w.rhs_les, fef.w.fun);
 332     fef.w.app = fef.w.fun;
 333     paraview.write();
 334     
 335     /*
 336     update_surface();
 337     update_scalar_solution();
 338     rhs_and_solve_SPDE();
 339     // update_velocity(); // does not exists anymore
 340     error_on_intSurface(); // Error on surface(t_n)
 341     next_timeStep();
 342     Pattern_helper helper {*this};
 343     helper.finalize_scalarPDE_rhs();
 344     helper.solve_scalarPDE();
 345     // helper.errors_on_numSurface();
 346     // helper.plot_paraview();
 347     */
 348   }
 349 }
 350 void Brusselator_scheme::final_action(){
 351   fef.u.write(io.dgf_handler, "./final_");
 352   fef.w.write(io.dgf_handler, "./final_");
 353   fef.surface.write(io.dgf_handler, "./final_"); 
 354 }
 355 
 356 // ------------------------------------------------------------
 357 // Brusselator_scheme private
 358 
 359 void Brusselator_scheme::update_surface(){
 360   io.identity.interpolate(fef.surface.exact);
 361 
 362   // Save old surface for the velocity
 363   fef.surface.app = fef.surface.fun; 
 364 }
 365 void Brusselator_scheme::update_scalar_solution(){
 366   exact.u.interpolate(fef.u.exact);
 367   exact.w.interpolate(fef.w.exact);
 368 }
 369 
 370 void Brusselator_scheme::error_on_intSurface(){
 371   const auto dT = fix_grid.time_provider().deltaT();
 372   // scalar_error
 373   io.u << dT << '\t'
 374        << norm.l2_err(fef.u.exact, fef.u.fun) << '\t'
 375        << norm.h1_err(fef.u.exact, fef.u.fun) << std::endl;
 376   io.w << dT << '\t'
 377        << norm.l2_err(fef.w.exact, fef.w.fun) << '\t'
 378        << norm.h1_err(fef.w.exact, fef.w.fun) << std::endl;
 379   // surface_error
 380   io.surface << dT << '\t'
 381              << norm.l2_err(fef.surface.exact, fef.surface.app)
 382              << '\t'
 383              << norm.h1_err(fef.surface.exact, fef.surface.app)
 384              << std::endl;
 385   // velocity error
 386   io.velocity << dT << '\t'
 387               << norm.l2_err(fef.velocity.fun, fef.velocity.exact)
 388               << '\t'
 389               << norm.h1_err(fef.velocity.fun, fef.velocity.exact)
 390               << std::endl;
 391 }
 392 
 393 void Brusselator_scheme::pre_loop_action(){
 394   PreLoop_helper helper {*this};
 395   helper.random_initialValues();
 396   // helper.analytic_initialValues();
 397   exact.X_ptr->interpolate(fef.surface.fun);
 398   helper.headLine_in_errFile();
 399   helper.save_surface();
 400   io.identity.interpolate(fef.surface.fun);
 401   // helper.plot_errors_in_errFile();
 402   // helper.plot_paraview();
 403   // next_timeStep();
 404   io.para << data << std::endl;
 405 }
 406 void Brusselator_scheme::rhs_and_solve_SPDE(){
 407   RhsAndSolve_helper helper {*this};
 408   helper.scalar_massMatrix();
 409   helper.brusselator_rhs();
 410   helper.addScaled_surfaceLoadVector();
 411   helper.solve_surface_and_save();
 412 }
 413 
 414 void Brusselator_scheme::test(){
 415   const auto it_end = pattern_timeSteps() + (data.last_step() >data.eps() ? 1 : 0 );
 416   std::cout << "t0: " << data.start_time() << '\n'
 417             << "dT: " << data.global_timeStep() << '\n'
 418             << "step no: " << it_end << '\n'
 419             << "last dT: " << data.last_step() << std::endl;
 420   std::cout << "time: " << fix_grid.time_provider().time() << std::endl;
 421   for(long it = 0; it < it_end; ++it){
 422     if(it < it_end-2){
 423       std::cout << "it < it_end-1" << std::endl;
 424       next_timeStep();
 425     }
 426     else{
 427       std::cout << "it >= it_end-1" << std::endl;
 428       fix_grid.next_timeStep(data.last_step());
 429     }
 430     std::cout << "time: " << fix_grid.time_provider().time() << std::endl;
 431   }
 432 }
 433 
 434 // ----------------------------------------------------------------------
 435 // Implementation of structs
 436 // Brusselator_scheme::Fef and Brusselator_scheme::Io
 437 
 438 Brusselator_scheme::Fef::Fef(const Esfem::Grid::Grid_and_time& gt)
 439   :u {"u", gt}, w {"w", gt},
 440    surface {"surface", gt}, velocity {"velocity", gt}
 441 {}
 442 
 443 Brusselator_scheme::Io::Io(const Esfem::Io::Parameter& p)
 444   :dgf_handler {p.grid()}, para {p}, u {"_u", p}, w {"_w", p},
 445    surface {"_X", p}, velocity {"_v",p}
 446 {}
 447 
 448 Brusselator_scheme::Init_data::Init_data(const Esfem::Grid::Grid_and_time& gt)
 449   :u {gt, Growth::promoting}, w {gt, Growth::inhibiting}, v {gt},
 450    X_ptr {// SecOrd_op::vIdata::new_ssef(gt)
 451      SecOrd_op::vIdata::new_sms(gt)}
 452 {}
 453 Brusselator_scheme::Init_data::Init_data(const Esfem::Grid::Grid_and_time& gt, const Esfem::Io::Parameter& p)
 454   :u {p, Growth::promoting}, w {p, Growth::inhibiting}, v {gt},
 455    X_ptr {SecOrd_op::vIdata::new_ssef(gt)}
 456 {}