1 /*! \file brusselator_algo_impl.h
2 \brief Implementation details for brusselator_algo.cpp
3
4 Revision history
5 --------------------------------------------------
6
7 Revised by Christian Power June 2016
8 Revised by Christian Power April 2016
9 Revised by Christian Power March 2016
10 Originally written by Christian Power
11 (power22c@gmail.com) February 2016
12
13 Idea
14 --------------------------------------------------
15
16 Providing many helper classes for the Brusselator_scheme.
17 Since the we solve two scalar equation I have provided several
18 struct to avoid word repetition.
19
20 \author Christian Power
21 \date 6. June 2016
22 \copyright Copyright (c) 2016 Christian Power. All rights reserved.
23 */
24
25 #ifndef BRUSSELATOR_ALGO_IMPL_H
26 #define BRUSSELATOR_ALGO_IMPL_H
27
28 #include <string>
29 #include "esfem.h"
30 #include "brusselator_algo.h"
31
32 namespace Esfem{
33 //! For the EOC experiments
34 struct Rhs{
35 //! Right-hand side for \f$ u \f$
36 SecOrd_op::Rhs u;
37 //! Right-hand side for \f$ w \f$
38 SecOrd_op::Rhs w;
39 //! Get time provider
40 explicit Rhs(const Grid::Grid_and_time&);
41 // Rhs(const Io::Parameter&, const Grid::Grid_and_time&);
42 };
43
44 struct Scalar_solver{
45 SecOrd_op::Brusselator u; /*!< PDE with solver for u */
46 SecOrd_op::Brusselator w; /*!< PDE with solver for w */
47 Scalar_solver(const Io::Parameter&,
48 const Grid::Grid_and_time&,
49 const Grid::Scal_FEfun_set& u_set,
50 const Grid::Scal_FEfun_set& w_set);
51 /*!< \brief Use this constructor if you need to solve a system. */
52 Scalar_solver(const Io::Parameter&,
53 const Grid::Grid_and_time&,
54 const Grid::Scal_tiny_FEfun_set& u_set,
55 const Grid::Scal_tiny_FEfun_set& w_set);
56 /*!< \brief Use this constructor if you only need a mass matrix.
57 Never use this constructor if you need to solve a system.
58 */
59 };
60
61 // struct Err_cal{
62 // Esfem::Io::L2H1_calculator u;
63 // /*!< \brief Calculates the error norm for u. */
64 // Esfem::Io::L2H1_calculator w;
65 // /*!< \brief Calculates the error norm for w. */
66 // Err_cal(const Grid::Grid_and_time&,
67 // const Grid::Scal_FEfun_set& u_set,
68 // const Grid::Scal_FEfun_set& w_set);
69 // };
70 // /*!< \brief Class that calculates errors in the \f$L^2\f$-
71 // and \f$H^1\f$-norm.
72 // */
73
74 class PreLoop_helper{
75 public:
76 //! Get access to
77 explicit PreLoop_helper(Brusselator_scheme&);
78 /*!< \brief Modifies private data members of `Brusselator_scheme` */
79 //! Initial values via an analytic expression
80 /*! \sa Rhs */
81 void analytic_initialValues();
82 //! Initial values via a random distribution
83 /*! \sa Rhs */
84 void random_initialValues();
85 //! Save the initial surface into a temporary file.
86 void save_surface();
87 //! First line in the error file
88 void headLine_in_errFile();
89 //! Paraview output
90 void plot_paraview();
91 void prepare_rhs();
92 /*!< \brief Applies the mass matrix on the
93 old surface on `u` and `w`.
94 \warning Overwrites the value of rhs.
95 */
96 private:
97 //! Contains numerical solution
98 Brusselator_scheme& bs;
99 //! Reference to `bs.exact`
100 const Brusselator_scheme::Init_data& init_data;
101 // Err_cal err_cal; /*!< \brief Contains two output files. */
102 Io::Paraview paraview;
103 /*!< \brief Has reference to `bs.fef.u.fun` and
104 `bs.fef.w.fun`
105 */
106 Scalar_solver solver; /*!< \brief Brusselator solver */
107 };
108 class PrePattern_helper{
109 public:
110 explicit PrePattern_helper(Brusselator_scheme&);
111 /*!< \brief Modifies private data members of `Brusselator_scheme` */
112 // void finalize_rhs();
113 void rhs();
114 /*!< \brief Adds to member `rhs_les` from
115 `Brusselator_scheme::fef.u` and
116 `Brusselator_scheme::fef.w`
117
118 The new value of
119 `rhs_les` from `Brusselator_scheme::fef.u` respectively
120 `Brusselator_scheme::fef.w` will be
121 \f{gather*}{
122 (M\nodalValue{u})^n + \tau \gamma a M^{n+1} \nodalValue{1}, \\
123 (M\nodalValue{w})^n + \tau \gamma b M^{n+1} \nodalValue{1}.
124 \f}
125 Note that the surface is not changing at this stage.
126 */
127 void solve_pde();
128 /*!< \brief New value for member `fun` from
129 `Brusselator_scheme::fef.u` and
130 `Brusselator_scheme::fef.w`
131
132 We assume that the `rhs_les` has been prepared with
133 finalize_rhs(). The vaule of `fun` will be overwritten
134 with the solution of the following system for \f$u\f$ and \f$w\f$:
135 \f{gather*}{
136 (1+ \tau \gamma) (M\nodalValue{u})^{n+1} + \tau (A \nodalValue{u})^{n+1}
137 - \tau \gamma M^{n+1}(\nodalValue{u}^n, \nodalValue{w}^n)
138 \nodalValue{u}^{n+1} = \nodalValue{y}_{1}, \\
139 (M\nodalValue{w})^{n+1} + \tau D_c (A \nodalValue{w})^{n+1}
140 + \tau \gamma M^{n+1}(\nodalValue{u}^{n+1}, \nodalValue{u}^{n+1})
141 \nodalValue{w}^{n+1} = \nodalValue{y}_{2},
142 \f}
143 where \f$ \nodalValue{y}_{1,2}\f$ the right-hand side of the system is.
144 */
145 void prepare_rhs();
146 /*!< \brief New value for member `rhs_les` from
147 `Brusselator_scheme::fef.u` and
148 `Brusselator_scheme::fef.w`
149
150 `rhs_les` will have the following value for \f$u\f$ respectively
151 \f$w\f$:
152 \f{equation*}{
153 (M\nodalValue{u})^n \quad \lor\quad (M \nodalValue{w})^n.
154 \f}
155 */
156 void plot_paraview();
157 /*!< \brief Prints out `fun` of `Brusselator_scheme::fef.u` and
158 `Brusselator_scheme::fef.w`.
159 */
160 private:
161 Brusselator_scheme::Io& io;
162 /*!< \brief Reference to `Brusselator_scheme::io` */
163 Grid::Scal_FEfun_set& u;
164 /*!< \brief Reference to `Brusselator_scheme::fef.u` */
165 Grid::Scal_FEfun_set& w;
166 /*!< \brief Reference to `Brusselator_scheme::fef.w` */
167 const Dune::Fem::TimeProviderBase& tp;
168 /*!< From `Brusselator_scheme::fix_grid` */
169 // Err_cal err_cal; /*!< \brief Contains two output files. */
170 Io::Paraview paraview;
171 /*!< \brief Has reference to member `fun` from `u` and `w` */
172 Scalar_solver solver; /*!< \brief Brusselator solver */
173 };
174 /*!< \brief Implementation details for
175 Brusselator_scheme::prePattern_loop().
176 */
177 class RhsAndSolve_helper{
178 public:
179 //! Get private members
180 explicit RhsAndSolve_helper(Brusselator_scheme&);
181 //! Assemble \f$ (M\nodalValue{u})^n \f$ and \f$ (M\nodalValue{w})^n \f$
182 void scalar_massMatrix();
183 /*! \brief \f$
184 (M_3^n + \alpha A_3^n) \nodalValue{X}^n
185 + \tau \delta M_3^n(\nodalValue{u}^n,
186 \nodalValue{\surfaceNormal}) \f$
187 */
188 void brusselator_rhs();
189 //! Adds \f$ \tau G^n\f$
190 void addScaled_surfaceLoadVector();
191 /*! \brief \f$
192 \parentheses[\big]{M_3^n + (\alpha + \varepsilon\tau) A_3^n} \nodalValue{X}^{n+1}
193 = \nodalValue{Y}
194 \f$
195 */
196 void solve_surface_and_save();
197 private:
198 //! Solver for `X`
199 using Vector_solver = Esfem::SecOrd_op::Solution_driven;
200
201 //! \copybrief PreLoop_helper::bs
202 Brusselator_scheme& bs;
203 //! Reference to `bs.fef`
204 Brusselator_scheme::Fef& fef;
205 //! Temporally grid, hence `const`
206 const Grid::Grid_and_time grid;
207 //! `fef.u.fun` and `fef.u.rhs_les` on `grid`
208 Grid::Scal_tiny_FEfun_set u;
209 //! `fef.w.fun` and `fef.w.rhs_les` on `grid`
210 Grid::Scal_tiny_FEfun_set w;
211 //! `bs_fef.surface.fun` and `bs_fef.surface.rhs_les` on `grid`
212 Grid::Vec_tiny_FEfun_set X;
213 //! Mass matrix for `u` and `w`
214 Scalar_solver ss;
215 //! Solver for `X`
216 Vector_solver vs;
217 //! Right-hand side for the surface PDE
218 SecOrd_op::Vec_rhs v_rhs;
219 };
220 /*!< \brief Implementation details for
221 Brusselator_scheme::rhs_and_solve_SPDE().
222 */
223 class Pattern_helper{
224 public:
225 explicit Pattern_helper(Brusselator_scheme&);
226 /*!< \copydoc PrePattern_helper() */
227 void finalize_scalarPDE_rhs();
228 /*!< \copydoc PrePattern_helper::finalize_rhs() */
229 void solve_scalarPDE();
230 /*!< \copydoc PrePattern_helper::solve_pde() */
231 //! Interpolate exact solutions scalar and vector valued PDE.
232 /*! Roughly I interpolate the ESFEM solution on the exact grid and then
233 assign nodal values to the current grid.
234 */
235 void update_exactSolutions();
236 // \f$L^2\f$- and \f$H^1\f$-errors on the numerical surface
237 // Not used anymore // void errors_on_numSurface();
238 void plot_paraview();
239 /*!< \copydoc PrePattern_helper::plot_paraview() */
240 private:
241 Brusselator_scheme& bs; /*!< \copydoc PreLoop_helper::bs */
242 const Grid::Grid_and_time grid;
243 /*!< \copydoc RhsAndSolve_helper::grid brief */
244 Grid::Scal_FEfun_set u; /*!< \brief `fef.u` on `grid` */
245 Grid::Scal_FEfun_set w; /*!< \brief `fef.w` on `grid` */
246 //! Norms on the calculated Grid
247 Io::L2H1_calculator norm;
248 Io::Paraview paraview;
249 /*!< \brief Has reference to member `u.fun` and `w.fun`. */
250 Scalar_solver solver; /*!< \brief Solver for `u` and `w` */
251 //! Provides assembly and add_scaled methods.
252 Rhs load_vector;
253 };
254 /*!< \brief Implementation details for
255 Brusselator_scheme::pattern_loop().
256 */
257
258 // ----------------------------------------------------------------------
259 // helper functions
260
261 void interpolate(const SecOrd_op::Init_data&, Grid::FEfun_set<Grid::Scal_FEfun>&);
262 /*!< \brief Gives initial values for the members `fun`, `app` and `exact` */
263 void head_line(Esfem::Io::Error_stream&);
264 /*!< \brief The head line is `timestep L2err H1err` with tab alignment. */
265 // void write_error_line(Io::Error_stream& file, const Dune::Fem::TimeProviderBase&,
266 // const Io::L2H1_calculator&);
267 /*!< \brief Prints time step, L2 and H1 error to `file` with proper
268 tab alignment.
269 */
270
271 // ----------------------------------------------------------------------
272 // Inline implementation
273
274 inline void PreLoop_helper::plot_paraview(){
275 paraview.write();
276 }
277 inline void PrePattern_helper::plot_paraview(){
278 paraview.write();
279 }
280 inline void Pattern_helper::plot_paraview(){
281 paraview.write();
282 }
283
284 } // namespace Esfem
285
286 #endif // BRUSSELATOR_ALGO_IMPL_H