1 /*! \file brusselator_algo.h
2 \brief Numerical experiment for the solution driven paper
3
4 Revision history
5 --------------------------------------------------
6
7 Revised by Christian Power June 2016
8 Revised by Christian Power May 2016
9 Revised by Christian Power April 2016
10 Revised by Christian Power March 2016
11 Originally written by Christian Power
12 (power22c@gmail.com) February 2016
13
14 Idea
15 --------------------------------------------------
16
17 This header provides model classes and operator classes to solve
18 a tumor growth model proposed by Elliott and Styles via the ESFEM.
19 Enter path to writable directory in the macro variable FEF_PATH.
20
21 \author Christian Power
22 \date 7. June 2016
23 \copyright Copyright (c) 2016 Christian Power. All rights reserved.
24 */
25
26 #ifndef BRUSSELATOR_ALGO_H
27 #define BRUSSELATOR_ALGO_H
28
29 #include <memory>
30 #include <string>
31 #include "esfem.h"
32
33 #ifndef FEF_PATH
34 #error Give full path to folder in FEF_PATH
35 #endif
36
37 namespace Esfem{
38 //! ESFEM algorithm. Only this should be invoked by main.
39 /*! \param argc `argc` from `main`
40 \param[in] argv `argv` from `main`
41 */
42 void brusselator_algo(int argc, char** argv);
43
44 //! Implementation of the Elliott and Styles full discretization of the tumor problem
45 /*!
46 Partial differential equation
47 ==================================================
48
49 Parameter
50 --------------------------------------------------
51
52 - \f$a,b, D_{c} \in \R\f$ for the scalar equation with
53 \f$\gamma \sim \vol(\surface)\f$.
54 - \f$\varepsilon, \delta, \alpha \in \R\f$ for the surface equation, where
55 \f$\varepsilon\f$ and \f$\alpha\f$ are small regularization parameter.
56
57 Smooth problem
58 --------------------------------------------------
59
60 Find \f$u\colon \surface \to \R\f$ (growth-promoting),
61 \f$w\colon \surface \to \R\f$ (growth-inhibiting) and
62 \f$X\colon \surface_{0} \times [0,T] \to \R^{m+1}\f$ such that
63 \f{gather*}{
64 \matd u + u \diver(v) - \laplaceBeltrami u = f_1(u,w) + f_{(1)}, \\
65 \matd w + w \diver(v) - D_c \laplaceBeltrami w = f_2(u,w) + f_{(2)},
66 \f}
67 where for \f$f_1,\, f_2\f$ we use the Brusselator model
68 \f{equation*}{
69 f_1(u,w) = \gamma (a - u + u^2 w) \quad \land \quad f_2(u,w)
70 = \gamma (b - u^2 w),
71 \f}
72 and \f$ f_{(1)}\f$ and \f$ f_{(2)}\f$ are some forcing terms,
73 and for the surface
74 \f{align*}{
75 v - \alpha \laplaceBeltrami v = {} &
76 \parentheses[\big]{\varepsilon (-\meanCurvature) + \delta u} \surfaceNormal
77 + g
78 = \varepsilon \laplaceBeltrami X + \delta u \surfaceNormal + g, \\
79 \dell_{t} X = {} & v(X),
80 \f}
81 where \f$ g\f$ is a forcing term.
82
83 Finite element discretization
84 --------------------------------------------------
85
86 Find \f$\nodalValue{u}\colon I \to \R^{N}\f$ (growth-promoting nodal values),
87 \f$\nodalValue{w}\colon I \to \R^{N}\f$ (growth-inhibiting
88 nodal values) and
89 \f$\nodalValue{X}\colon I \to \R^{3N}\f$ (surface nodal values) such that
90 \f{gather*}{
91 \parentheses[\big]{M(X) + \alpha A(X)} \dell_t X =
92 -\varepsilon A(X)X + \delta M(u,\nodalValue{\surfaceNormal})
93 + G, \\
94 \dell_t \parentheses[\big]{M(X) \nodalValue{u} } + A(X) \nodalValue{u}
95 = \gamma \parentheses[\big]{a M(X) \nodalValue{1}
96 - M(X)\nodalValue{u}
97 + M(X; \nodalValue{u},\nodalValue{w}) \nodalValue{u}}
98 + F_{(1)} \\
99 \dell_t \parentheses[\big]{M(X) \nodalValue{w}}
100 + D_c A(X) \nodalValue{w}
101 = \gamma \parentheses[\big]{b M(X) \nodalValue{1}
102 - M(X; \nodalValue{u},\nodalValue{u}) \nodalValue{w}}
103 + F_{(2)}.
104 \f}
105 We note that instead of the \f$ L^2\f$-projection we use the interpolation
106 of \f$ g\f$, \f$ f_{(1)}\f$ respectively \f$ f_{(2)}\f$.
107
108 Full discretization (Elliott+Styles discretization)
109 --------------------------------------------------
110
111 We perform three steps.
112
113 1. Given \f$\nodalValue{X}^n,\, \nodalValue{u}^n,\, \nodalValue{w}^n\f$
114 solve for \f$\nodalValue{X}^{n+1}\f$
115 \f{equation*}{
116 \parentheses[\big]{M_3^n + (\alpha + \varepsilon\tau) A_3^n}
117 \nodalValue{X}^{n+1}
118 = (M_3^n + \alpha A_3^n) \nodalValue{X}^n
119 + \tau \delta M_3^n(\nodalValue{u}^n,
120 \nodalValue{\surfaceNormal}) + \tau G^n,
121 \f}
122 where \f$\nodalValue{\surfaceNormal}^n\f$ is elementwise normal.
123 2. Given \f$\nodalValue{X}^{n+1},\, \nodalValue{u}^n,\, \nodalValue{w}^n\f$
124 solve for \f$\nodalValue{u}^{n+1}\f$
125 \f{equation*}{
126 (1+ \tau \gamma) (M\nodalValue{u})^{n+1} + \tau (A \nodalValue{u})^{n+1}
127 - \tau \gamma M^{n+1}(\nodalValue{u}^n, \nodalValue{w}^n)
128 \nodalValue{u}^{n+1}
129 = (M\nodalValue{u})^n + \tau \gamma a M^{n+1} \nodalValue{1}
130 + \tau F^{n+1}_{(1)},
131 \f}
132 where \f$M(a,b)\f$ a \f$4\f$ tensor is, namely
133 \f$M_{ijkl} = \int \chi_i \chi_j \chi_k \chi_l\f$
134 and \f$\nodalValue{1}\f$ means the
135 constant \f$1\f$ finite element function.
136 3. Given \f$\nodalValue{X}^{n+1},\, \nodalValue{u}^{n+1},\,
137 \nodalValue{w}^n\f$
138 solve for \f$\nodalValue{w}^{n+1}\f$
139 \f{equation*}{
140 (M\nodalValue{w})^{n+1} + \tau D_c (A \nodalValue{w})^{n+1}
141 + \tau \gamma M^{n+1}(\nodalValue{u}^{n+1}, \nodalValue{u}^{n+1})
142 \nodalValue{w}^{n+1}
143 = (M \nodalValue{w})^n + \tau \gamma b M^{n+1} \nodalValue{1}
144 + \tau F^{n+1}_{(2)}.
145 \f}
146 */
147 class Brusselator_scheme{
148 public:
149 Brusselator_scheme(int argc, char** argv,
150 const std::string& parameter_fname);
151 /*!< \brief The constructor that also performs the
152 first part before the loop enters
153 \param argc `argc` from `main`
154 \param[in] argv `argv` from `main`
155 \param parameter_fname Preferable absolute path to parameter file.
156 \warning Absolute path differs on different operating systems.
157 */
158
159 /*! \name Prescribed ESFEM */
160 //@{
161 //! Standard Dziuk Elliott evolving surface finite element method
162 void standard_esfem();
163 //@}
164
165 //! Experiment, where right-hand side is calculated for a known solution
166 /*! As exact solution for the scalar valued surface equation I chose
167 \f{equation*}{
168 u(x,y,z,t) = x y e^{-6t} \quad \text{and}\quad
169 w(x,y,z,t) = y z e^{-6t}.
170 \f}
171 For the surface evolution I chose
172 \f{equation*}{
173 \Phi(x,t) := r(t) x, \quad
174 r(t) := \frac{r_{end} r_0}{r_{end} e^{-kt} + r_0 (1-e^{-kt})},
175 \f}
176 where \f$ r(t)\f$ is the logistic growth function, which satisfies the
177 following ODE
178 \f{equation*}{
179 \dot{r} = k \Bigl( 1 - \frac{r}{r_{end}}\Bigr) r,\quad r(0) = r_0.
180 \f}
181 From this it follows easily that the velocity is given by
182 \f{equation*}{
183 v(x,t) = k \Bigl(1 - \frac{r}{r_{end}}\Bigr) x.
184 \f}
185 */
186 void eoc_logisticSphere();
187
188 //! Dziuk mean curvature flow ESFEM
189 /*! (This is not any more true) Experiment reads as follows:
190 Stationary surface \f$\surface = S^2\f$
191 with exact solution \f$X=e^{-t}(xy, yz, xz)\f$, with the formula
192 \f[
193 \Delta f = \Delta_{\R^{n+1}} f - H D_{\R^{n+1}} f \cdot n
194 - D^2_{\R^{n+1}} f(n,n),
195 \f]
196 and
197 \f[
198 H = div(n) = \frac{n}{|x|}
199 \f]
200 one easily sees
201 \f[
202 \Delta(xy) = -\frac{2n+2}{|x|^2} xy,
203 \f]
204 which implies that the right-hand side of the surface PDE vanishes.
205 */
206 void eoc_mcf();
207
208 //! Sphere Dalquist
209 /*! Initial surface is unit sphere in \f$\R^3\f$. Exact solution is
210 \f[
211 \Phi(x,t) := r(t)x,\quad r(t):= e^t
212 \f]
213 The PDE is
214 \f[
215 v - \alpha \Delta v - \varepsilon \Delta x
216 = f(x,t) = (1 + (\alpha + \varepsilon) 2 e^{-2t})x.
217 \f]
218 */
219 void sd();
220 //! Surface logistic sphere
221 /*! As exact solution for the surface evolution I choose
222 \f{equation*}{
223 \Phi(x,t) := r(t) x, \quad
224 r(t) := \frac{r_{end} r_0}{r_{end} e^{-kt} + r_0 (1-e^{-kt})}.
225 \f}
226 \f$r(t)\f$ satiesfies the ODE
227 \f[
228 \dot{r}(t) = k \left( 1 - \frac{r(t)}{r_{end}}\right) r(t),
229 \f]
230 which implies for the velocity
231 \f[
232 v(x,t) = k \left( 1 - \frac{r(t)}{r_{end}} \right)x
233 =: \tilde{a} x
234 \f]
235 I do not consider coupling. For the surface PDE I choose
236 \f[
237 v - \alpha \Delta v - \varepsilon \Delta x - \delta u \vec{n} = f,
238 \f]
239 where \f$f\f$ must be
240 \f[
241 f = \left( \tilde{a} +
242 \frac{(\alpha \tilde{a} + \varepsilon)H - \delta u}{|x|}\right) x,
243 \f]
244 where we have used \f$ \Delta x = - H \vec{n}\f$, where \f$H\f$
245 is the mean curvature
246 (without aritmetic mean) and \f$\vec{n}\f$ is the outwards pointing normal, and
247 \f$ \vec{n} = \frac{x}{|x|}\f$. It holds \f$H = \frac{n}{|x|}\f$,
248 where \f$n\f$ is the dimension of the sphere. Note that \f$|x| = r(t)\f$
249 on the exact surface.
250 \sa Esfem::SecOrd_op::vRhs::new_sls(), Esfem::SecOrd_op::vIdata::new_sls()
251 */
252 void eoc_sls();
253
254 //! EOC experiment for solution driven paper
255 /*! As exact solution for the surface evolution I choose
256 \f{equation*}{
257 \Phi(x,t) := r(t) x, \quad
258 r(t) := \frac{r_{end} r_0}{r_{end} e^{-kt} + r_0 (1-e^{-kt})}.
259 \f}
260 \f$r(t)\f$ satiesfies the ODE
261 \f[
262 \dot{r}(t) = k \left( 1 - \frac{r(t)}{r_{end}}\right) r(t),
263 \f]
264 which implies for the velocity
265 \f[
266 v(x,t) = k \left( 1 - \frac{r(t)}{r_{end}} \right)x
267 =: \tilde{a} x
268 \f]
269 As exact solution for the scalar diffusion equation I choose
270 \f[
271 u(x,t) := x_1 x_2 e^{-6t}.
272 \f]
273 The coupled PDE equations reads as
274 \f{gather*}{
275 \partial^{\bullet} u + u \nabla\cdot v - \Delta u = f \\
276 v - \Delta v - \Delta X - \delta u = g,
277 \f}
278 where \f$g\f$ is given via
279 \f[
280 g = \left( \tilde{a} +
281 \frac{(\alpha \tilde{a} + \varepsilon)H - \delta u}{|x|}\right) x,
282 \f]
283 and the complicated \f$f\f$ is computed via Sage.
284 \pre Parameter \f$\gamma\f$ has to be zero.
285 \sa eoc_sls(), Esfem::SecOrd_op::vRhs::new_sls(),
286 Esfem::SecOrd_op::sRhs::new_sdp_u(), Esfem::SecOrd_op::vIdata::new_sls(),
287 Esfem::SecOrd_op::sIdata::new_1ssef()
288 */
289 void eoc_sdp();
290
291 //! Run a simple test
292 void test();
293
294 /*! \name Loop action */
295 //@{
296 void prePattern_loop();
297 /*!< \brief In some sense calculates the inital data for
298 the solution driven problem. Starts from \f$t_0\f$.
299 */
300 void intermediate_action();
301 /*!< \brief To be used between prePattern_loop()
302 and pattern_loop().
303
304 The inital data has been created. It will be saved
305 and the right hand side for the surface PDE will be created.
306 \warning Do not use next_timeStep() afterwards.
307 */
308 void pattern_loop();
309 /*!< \brief To be used in the second for-loop.
310 If prePattern_loop() is off, then \f$t_0\f$ is here.
311
312 At this stage the tumor is growing.
313 \warning Do not use next_timeStep() afterwards.
314 */
315 void final_action();
316 /*!< \brief To be used after the second for-loop to save some data. */
317 //@}
318 private:
319 struct Io{
320 //! Functor for saving the current grid.
321 SecOrd_op::Identity identity {};
322
323 const Esfem::Io::Dgf::Handler dgf_handler;
324 /*!< \brief Converts finite element function into dgf file. */
325
326 //! File to capture PDE parameter
327 Esfem::Io::Error_stream para;
328 Esfem::Io::Error_stream u;
329 /*!< \brief File to record errors of u. */
330 Esfem::Io::Error_stream w;
331 /*!< \brief File to record errors of w. */
332 //! File to plot errors in the surface
333 Esfem::Io::Error_stream surface;
334 //! File to plot errors in the velocity
335 Esfem::Io::Error_stream velocity;
336
337 //! Get file name for error streams.
338 Io(const Esfem::Io::Parameter&);
339 };
340 /*!< \brief Shortens Brusselator_scheme().
341
342 Members are used for input and output of
343 the nodal values of the finite element functions
344 in `fef` or to calculate \f$L^2\f$- or \f$H^1\f$-norms of the errors.
345 */
346 //! Initial data respectively exact solution for numerical experiments
347 struct Init_data{
348 //! Initial data respectively exact solution functor for \f$u\f$
349 SecOrd_op::Init_data u;
350 //! Initial data respectively exact solution functor for \f$w\f$
351 SecOrd_op::Init_data w;
352 //! Interpolation functor for the exact velocity
353 SecOrd_op::Exact_velocity v;
354 //! Exact solution of the surface
355 std::unique_ptr<SecOrd_op::vIdata> X_ptr;
356
357 //! Provides analytically given initial data.
358 explicit Init_data(const Grid::Grid_and_time&);
359 //! Provides uniform distributed random inital values.
360 Init_data(const Grid::Grid_and_time&, const Esfem::Io::Parameter&);
361 };
362 struct Fef{
363 Grid::Scal_FEfun_set u;
364 /*!< \brief Container for growth promoting numerical solution */
365 Grid::Scal_FEfun_set w;
366 /*!< \brief Container for growth inhibiting numerical solution */
367 Grid::Vec_FEfun_set surface;
368 /*!< \brief Container for the numerical solution of the surface
369 \warning Only this container is allowed to write into and read from
370 a dgf file.
371 */
372 //! Analytically given exact velocity.
373 Grid::Vec_FEfun_set velocity;
374 const std::string tmpFile_path {FEF_PATH};
375 /*!< \brief Directory which I have read and write access.
376 \warning `FEF_PATH` is a macro variable which has be set by
377 the makefile.
378 */
379 Fef(const Grid::Grid_and_time&);
380 };
381 /*!< \brief Shortens Brusselator_scheme().
382
383 Collects all finite element functions and serves as a
384 backup container.
385 \todo Add error checking for `tmpFile_path`.
386 */
387
388 Esfem::Io::Parameter data;
389 /*!< \brief Contains parameter from `tumor_parameter.txt`. */
390 //! Error streams and identity functor for the surface
391 Io io;
392 //! Analytically given grid with absolute time provider.
393 Esfem::Grid::Grid_and_time fix_grid;
394 //! Norms on the analytically given grid
395 Esfem::Io::L2H1_calculator norm;
396 //! Finite element functions
397 Fef fef;
398 //! Exact solution for \f$u\f$, \f$w\f$ and \f$v\f$
399 Init_data exact;
400
401 // ------------------------------------------------------------
402 // Helper member functions
403
404 /*! \name Helper classes for the for-loops */
405 //@{
406 friend class PreLoop_helper;
407 /*!< \brief Used in pre_loop_action(). */
408 friend class PrePattern_helper;
409 /*!< \brief Used in prePattern_loop(). */
410 friend class Pattern_helper;
411 /*!< \brief Used in pattern_loop(). */
412 friend class RhsAndSolve_helper;
413 /*!< \brief Used in rhs_and_solve_SPDE(). */
414 //@}
415
416 //! Assign a new value to `fef.surface.exact` and `fef.surface.app`
417 void update_surface();
418 //! Assign a new value to `fef.u.exact` and `fef.w.exact`
419 void update_scalar_solution();
420
421 //! Lifted printing
422 template<class F>
423 void print(Esfem::Io::Error_stream& os, const F& fem){
424 os << fix_grid.time_provider().deltaT() << ' '
425 << norm.l2_err(fem.fun, fem.exact) << ' '
426 << norm.h1_err(fem.fun, fem.exact) << std::endl;
427 }
428
429 //! Calculate the velocity via simple differential quotient
430 /*! \param[in] xn_first Iterator to first point of the new surface
431 \param[in] xn_last Iterator to last point of the new surface
432 \param[in] xo_first Iterator to first point of the previous surface
433 \param[out] v_first Iterator to the first value of the velocity
434 \pre Both `xo_first` and `v_first` point to as many elements as
435 `xn_first` does.
436 */
437 template<class It1, class It2>
438 void calculate_velocity(It1 xn_first, It1 xn_last, It1 xo_first, It2 v_first){
439 const double dT = fix_grid.time_provider().deltaT();
440 for(; xn_first != xn_last; ++xn_first, ++xo_first, ++v_first)
441 *v_first = (*xn_first - *xo_first)/dT;
442 }
443
444
445 //! Plot error of `fef` on the interpolated surface
446 /*! \pre The exact solution has been updated.
447 \sa update_exact_surface(), update_exact_velocity(),
448 update_scalar_solution()
449 */
450 void error_on_intSurface();
451
452 //! Constructor helper
453 /*! \pre Should only be invoked by Brusselator_scheme().*/
454 void pre_loop_action();
455 //! Helper for pattern_loop().
456 /*! Generates first part of the right-hand side for the scalar
457 SPDE. Then solves the vector SPDE and prints out a dgf file.
458 */
459 void rhs_and_solve_SPDE();
460
461 /*! \name Flow control */
462 //@{
463 void next_timeStep();
464 /*!< \brief Increments the next time step in `fix_grid`. */
465 long prePattern_timeSteps() const;
466 /*!< \brief Maximum number of time steps for prePattern_loop(). */
467 long pattern_timeSteps() const;
468 /*!< \brief Maximum number of time steps for pattern_loop(). */
469 //! Absolute number of time steps
470 long time_steps() const;
471 //@}
472 };
473
474 // ----------------------------------------------------------------------
475 // Inline implementation
476
477 inline void Brusselator_scheme::next_timeStep(){
478 fix_grid.next_timeStep(data.global_timeStep());
479 }
480 inline long Brusselator_scheme::prePattern_timeSteps() const{
481 return data.prePattern_timeSteps();
482 }
483 inline long Brusselator_scheme::pattern_timeSteps() const{
484 return data.pattern_timeSteps();
485 }
486 inline long Brusselator_scheme::time_steps() const{
487 return data.max_timeSteps();
488 }
489 }
490
491 #endif // BRUSSELATOR_ALGO_H