root/src/brusselator_algo_impl.cpp

DEFINITIONS

1. analytic_initialValues
2. random_initialValues
3. save_surface
4. headLine_in_errFile
5. prepare_rhs
6. rhs
7. solve_pde
8. prepare_rhs
9. scalar_massMatrix
10. brusselator_rhs
11. addScaled_surfaceLoadVector
12. solve_surface_and_save
13. finalize_scalarPDE_rhs
14. solve_scalarPDE
15. update_exactSolutions
16. interpolate
17. head_line

   1 /*! \file brusselator_algo_impl.cpp
   2     \brief Implementation details for brusselator_algo_impl.h
   3 
   4      Revision history
   5      --------------------------------------------------
   6 
   7           Revised by Christian Power April 2016
   8           Revised by Christian Power March 2016
   9           Originally written by Christian Power
  10                (power22c@gmail.com) February 2016
  11 
  12      Idea
  13      --------------------------------------------------
  14 
  15      Providing many struct, such that functions like initial value,
  16      right-hand side function or PDE solver do not have to be specified for 
  17      \f$ u\f$ and \f$ w\f$ separately.  
  18 
  19     \author Christian Power
  20     \date 25. April 2016
  21     \copyright Copyright (c) 2016 Christian Power.  All rights reserved.
  22 */
  23 
  24 #include "brusselator_algo_impl.h"
  25 #include "esfem_error.h"
  26 
  27 //! Reporting errors in this file with class.
  28 using Bruss_error = Esfem::BrusselatorScheme_error;
  29 //! Four finite element functions of type \f$\surface_h \to \R\f$
  30 using Scal_FEfun_set = Esfem::Grid::Scal_FEfun_set;
  31 //! Four finite element functions of type \f$\surface_h \to \R^3\f$
  32 using Vec_FEfun_set = Esfem::Grid::Vec_FEfun_set;
  33 
  34 using Esfem::SecOrd_op::Identity;
  35 using Esfem::Rhs;
  36 using Esfem::Scalar_solver;
  37 using Esfem::PrePattern_helper;
  38 using Esfem::PreLoop_helper;
  39 using Esfem::RhsAndSolve_helper;
  40 using Esfem::Pattern_helper;
  41 
  42 // ----------------------------------------------------------------------
  43 // Implementation of structs 
  44 
  45 Rhs::Rhs(const Grid::Grid_and_time& gt)
  46   :u {gt, Growth::promoting}, w {gt, Growth::inhibiting}
  47 {}
  48 
  49 Scalar_solver::Scalar_solver
  50 (const Esfem::Io::Parameter& p,
  51  const Esfem::Grid::Grid_and_time& g,
  52  const Esfem::Grid::Scal_FEfun_set& u_set,
  53  const Esfem::Grid::Scal_FEfun_set& w_set)
  54   :u {p, g, Growth::promoting, u_set.app, w_set.app},
  55    w {p, g, Growth::inhibiting, u_set.app, u_set.app}
  56 {}
  57 
  58 Scalar_solver::Scalar_solver
  59 (const Esfem::Io::Parameter& p,
  60  const Esfem::Grid::Grid_and_time& g,
  61  const Esfem::Grid::Scal_tiny_FEfun_set& u_set,
  62  const Esfem::Grid::Scal_tiny_FEfun_set& w_set)
  63   :u {p, g, Growth::promoting, u_set.fun, w_set.fun},
  64    w {p, g, Growth::inhibiting, u_set.fun, u_set.fun}
  65 {}
  66 
  67 // Err_cal::Err_cal(const Esfem::Grid::Grid_and_time& g,
  68 //               const Scal_FEfun_set& u_set,
  69 //               const Scal_FEfun_set& w_set)
  70 //   :u {g, u_set.exact, u_set.fun},
  71 //   w {g, w_set.exact, w_set.fun}
  72 // {}
  73 
  74 // ----------------------------------------------------------------------
  75 // Implementation PreLoop_helper
  76 
  77 PreLoop_helper::PreLoop_helper(Brusselator_scheme& bs_input)
  78   :bs {bs_input},
  79    init_data {bs.exact},
  80    paraview {bs.data, bs.fix_grid, bs.fef.u.fun, bs.fef.w.fun},
  81    solver {bs.data, bs.fix_grid, bs.fef.u, bs.fef.w}
  82 {}
  83 
  84 void PreLoop_helper::analytic_initialValues(){
  85   interpolate(bs.exact.u, bs.fef.u);
  86   interpolate(bs.exact.w, bs.fef.w);
  87   // bs.io.identity.interpolate(bs.fef.surface.fun);
  88 }
  89 void PreLoop_helper::random_initialValues(){
  90   interpolate(init_data.u, bs.fef.u);
  91   interpolate(init_data.w, bs.fef.w);
  92   bs.io.identity.interpolate(bs.fef.surface.fun);
  93 }
  94 void PreLoop_helper::save_surface(){
  95   bs.fef.surface.write(bs.io.dgf_handler, bs.fef.tmpFile_path);
  96 }
  97 void PreLoop_helper::headLine_in_errFile(){
  98   head_line(bs.io.u);
  99   head_line(bs.io.w);
 100   head_line(bs.io.surface);
 101   head_line(bs.io.velocity);
 102 }
 103 void PreLoop_helper::prepare_rhs(){
 104   auto& u = bs.fef.u;
 105   auto& w = bs.fef.w;
 106   solver.u.mass_matrix(u.fun, u.rhs_les);
 107   solver.w.mass_matrix(w.fun, w.rhs_les);
 108 }
 109 
 110 // ----------------------------------------------------------------------
 111 // Implementation PrePattern_helper 
 112 
 113 PrePattern_helper::PrePattern_helper(Brusselator_scheme& bs)
 114   :io {bs.io},
 115   u {bs.fef.u},
 116   w {bs.fef.w},
 117   tp {bs.fix_grid.time_provider()},
 118   paraview {bs.data, bs.fix_grid, u.fun, w.fun},
 119   solver {bs.data, bs.fix_grid, u, w}
 120 {}
 121 
 122 // void PrePattern_helper::finalize_rhs(){
 123 void PrePattern_helper::rhs(){
 124   solver.u.mass_matrix(u.fun, u.rhs_les);
 125   solver.w.mass_matrix(w.fun, w.rhs_les);
 126   solver.u.add_massMatrixConstOne_to(u.rhs_les);
 127   solver.w.add_massMatrixConstOne_to(w.rhs_les);
 128 }
 129 void PrePattern_helper::solve_pde(){
 130   solver.u.solve(u.rhs_les, u.fun);
 131   u.app = u.fun;    
 132   solver.w.solve(w.rhs_les, w.fun);
 133   w.app = w.fun;    
 134 }
 135 void PrePattern_helper::prepare_rhs(){
 136   solver.u.mass_matrix(u.fun, u.rhs_les);
 137   solver.w.mass_matrix(w.fun, w.rhs_les);
 138 }
 139 
 140 // ----------------------------------------------------------------------
 141 // Implementation RhsAndSolve_helper
 142 
 143 RhsAndSolve_helper::RhsAndSolve_helper(Brusselator_scheme& bs_input)
 144 try :bs {bs_input},
 145   fef {bs.fef},
 146   grid {bs.data,
 147       Grid::compose_dgfName(fef.surface.fun.name(), fef.tmpFile_path), 
 148       bs.fix_grid.time_provider().time()},
 149   u {fef.u, grid},
 150   w {fef.w, grid},
 151   X {fef.surface, grid},
 152   ss {bs.data, grid, u, w},
 153   vs {bs.data, grid, u.fun},
 154   v_rhs {bs.fix_grid}
 155 {}
 156 catch(...){
 157   std::throw_with_nested(Bruss_error {"RhsAndSolve_helper()"});
 158  }
 159 void RhsAndSolve_helper::scalar_massMatrix(){
 160   ss.u.mass_matrix(u.fun, u.rhs_les);
 161   ss.w.mass_matrix(w.fun, w.rhs_les);
 162   fef.u.rhs_les = u.rhs_les;
 163   fef.w.rhs_les = w.rhs_les;
 164 }
 165 void RhsAndSolve_helper::brusselator_rhs(){
 166   vs.brusselator_rhs(X.fun, X.rhs_les); 
 167 }
 168 void RhsAndSolve_helper::addScaled_surfaceLoadVector(){
 169   v_rhs.assemble_and_addScaled_to(X.rhs_les);
 170 }
 171 void RhsAndSolve_helper::solve_surface_and_save(){
 172   vs.solve(X.rhs_les, X.fun);
 173   fef.surface.fun = X.fun;
 174   fef.surface.write(bs.io.dgf_handler, fef.tmpFile_path);
 175 }
 176 
 177 // ----------------------------------------------------------------------
 178 // Implementation Pattern_helper
 179 
 180 /*! \todo scalar valued rhs is missing */
 181 Pattern_helper::Pattern_helper(Brusselator_scheme& bs_input)
 182   :bs {bs_input},
 183    grid
 184   {bs.data,
 185       Grid::compose_dgfName(bs.fef.surface.fun.name(), bs.fef.tmpFile_path ), 
 186       bs.fix_grid.time_provider().time()},
 187   u {bs.fef.u, grid},
 188   w {bs.fef.w, grid},
 189   norm {grid},
 190   paraview {bs.data, grid, u.fun, w.fun},
 191   solver {bs.data, grid, u, w},
 192   load_vector {grid}
 193 {}
 194 
 195 void Pattern_helper::finalize_scalarPDE_rhs(){
 196   solver.u.add_massMatrixConstOne_to(u.rhs_les);
 197   load_vector.u.assemble_and_addScaled_to(u.rhs_les);
 198   solver.w.add_massMatrixConstOne_to(w.rhs_les);
 199   load_vector.w.assemble_and_addScaled_to(w.rhs_les);
 200 }
 201 void Pattern_helper::solve_scalarPDE(){
 202   solver.u.solve(u.rhs_les, u.fun);
 203   u.app = u.fun;    
 204   solver.w.solve(w.rhs_les, w.fun);
 205   w.app = w.fun;
 206   bs.fef.u = u;
 207   bs.fef.w = w;
 208 }
 209 void Pattern_helper::update_exactSolutions(){
 210   // bs.exact.u.interpolate(bs.fef.u.exact);
 211   // bs.exact.w.interpolate(bs.fef.w.exact);
 212   bs.exact.u.interpolate(u.exact);
 213   bs.exact.w.interpolate(w.exact);
 214 }
 215 // void Pattern_helper::errors_on_numSurface(){
 216 //   const auto& tp = grid.time_provider();
 217 //   write_error_line(bs.io.u, tp, err_cal.u);
 218 //   write_error_line(bs.io.w, tp, err_cal.w);  
 219 // }
 220 
 221 // ----------------------------------------------------------------------
 222 // helper functions
 223 
 224 void Esfem::interpolate(const SecOrd_op::Init_data& id, Scal_FEfun_set& f){
 225   id.interpolate(f.fun);
 226   f.app = f.fun;
 227   f.exact = f.fun;
 228 }
 229 void Esfem::head_line(Io::Error_stream& file){
 230   file << "timestep" << "\t"
 231        << "L2err" << "\t\t"
 232        << "H1err" << std::endl;
 233   file << std::scientific;
 234 }
 235 // void Esfem::write_error_line(Io::Error_stream& file,
 236 //                           const Dune::Fem::TimeProviderBase& tp,
 237 //                           const Io::L2H1_calculator& cal){
 238 //   file << tp.deltaT() << '\t'
 239 //        << cal.l2_err() << '\t'
 240 //        << cal.h1_err() << std::endl; 
 241 // }