Implementation of the Elliott and Styles full discretization of the tumor problem. More...

#include <brusselator_algo.h>

Collaboration diagram for Esfem::Brusselator_scheme:

[legend]

Classes
struct	Fef
	Shortens Brusselator_scheme(). More...

struct	Init_data
	Initial data respectively exact solution for numerical experiments. More...

struct	Io
	Shortens Brusselator_scheme(). More...

Public Member Functions
	Brusselator_scheme (int argc, char **argv, const std::string &parameter_fname)
	The constructor that also performs the first part before the loop enters. More...

void	eoc_logisticSphere ()
	Experiment, where right-hand side is calculated for a known solution. More...

void	eoc_mcf ()
	Dziuk mean curvature flow ESFEM. More...

void	sd ()
	Sphere Dalquist. More...

void	eoc_sls ()
	Surface logistic sphere. More...

void	eoc_sdp ()
	EOC experiment for solution driven paper. More...

void	test ()
	Run a simple test.

Prescribed ESFEM
void	standard_esfem ()
	Standard Dziuk Elliott evolving surface finite element method.

Loop action
void	prePattern_loop ()
	In some sense calculates the inital data for the solution driven problem. Starts from .

void	intermediate_action ()
	To be used between prePattern_loop() and pattern_loop(). More...

void	pattern_loop ()
	To be used in the second for-loop. If prePattern_loop() is off, then is here. More...

void	final_action ()
	To be used after the second for-loop to save some data.

Private Member Functions
void	update_surface ()
	Assign a new value to `fef.surface.exact` and `fef.surface.app`

void	update_scalar_solution ()
	Assign a new value to `fef.u.exact` and `fef.w.exact`

template<class F >
void	print (Esfem::Io::Error_stream &os, const F &fem)
	Lifted printing.

template<class It1 , class It2 >
void	calculate_velocity (It1 xn_first, It1 xn_last, It1 xo_first, It2 v_first)
	Calculate the velocity via simple differential quotient. More...

void	error_on_intSurface ()
	Plot error of `fef` on the interpolated surface. More...

void	pre_loop_action ()
	Constructor helper. More...

void	rhs_and_solve_SPDE ()
	Helper for pattern_loop(). More...

Flow control
void	next_timeStep ()
	Increments the next time step in `fix_grid`.

long	prePattern_timeSteps () const
	Maximum number of time steps for prePattern_loop().

long	pattern_timeSteps () const
	Maximum number of time steps for pattern_loop().

long	time_steps () const
	Absolute number of time steps.

Private Attributes
Esfem::Io::Parameter	data
	Contains parameter from `tumor_parameter.txt`.

Io	io
	Error streams and identity functor for the surface.

Esfem::Grid::Grid_and_time	fix_grid
	Analytically given grid with absolute time provider.

Esfem::Io::L2H1_calculator	norm
	Norms on the analytically given grid.

Fef	fef
	Finite element functions.

Init_data	exact
	Exact solution for , and .

Friends
Helper classes for the for-loops
class	PreLoop_helper
	Used in pre_loop_action().

class	PrePattern_helper
	Used in prePattern_loop().

class	Pattern_helper
	Used in pattern_loop().

class	RhsAndSolve_helper
	Used in rhs_and_solve_SPDE().

Detailed Description

Implementation of the Elliott and Styles full discretization of the tumor problem.

Partial differential equation

Parameter

$a,b, D_{c} \in \R$ for the scalar equation with $\gamma \sim \vol(\surface)$ .
$\varepsilon, \delta, \alpha \in \R$ for the surface equation, where $\varepsilon$ and $\alpha$ are small regularization parameter.

Smooth problem

Find $u\colon \surface \to \R$ (growth-promoting), $w\colon \surface \to \R$ (growth-inhibiting) and $X\colon \surface_{0} \times [0,T] \to \R^{m+1}$ such that

$\begin{gather*} \matd u + u \diver(v) - \laplaceBeltrami u = f_1(u,w) + f_{(1)}, \\ \matd w + w \diver(v) - D_c \laplaceBeltrami w = f_2(u,w) + f_{(2)}, \end{gather*}$

where for $f_1,\, f_2$ we use the Brusselator model

$\begin{equation*} f_1(u,w) = \gamma (a - u + u^2 w) \quad \land \quad f_2(u,w) = \gamma (b - u^2 w), \end{equation*}$

and $f_{(1)}$ and $f_{(2)}$ are some forcing terms, and for the surface

$\begin{align*} v - \alpha \laplaceBeltrami v = {} & \parentheses[\big]{\varepsilon (-\meanCurvature) + \delta u} \surfaceNormal + g = \varepsilon \laplaceBeltrami X + \delta u \surfaceNormal + g, \\ \dell_{t} X = {} & v(X), \end{align*}$

where $ g$ is a forcing term.

Finite element discretization

Find $\nodalValue{u}\colon I \to \R^{N}$ (growth-promoting nodal values), $\nodalValue{w}\colon I \to \R^{N}$ (growth-inhibiting nodal values) and $\nodalValue{X}\colon I \to \R^{3N}$ (surface nodal values) such that

$\begin{gather*} \parentheses[\big]{M(X) + \alpha A(X)} \dell_t X = -\varepsilon A(X)X + \delta M(u,\nodalValue{\surfaceNormal}) + G, \\ \dell_t \parentheses[\big]{M(X) \nodalValue{u} } + A(X) \nodalValue{u} = \gamma \parentheses[\big]{a M(X) \nodalValue{1} - M(X)\nodalValue{u} + M(X; \nodalValue{u},\nodalValue{w}) \nodalValue{u}} + F_{(1)} \\ \dell_t \parentheses[\big]{M(X) \nodalValue{w}} + D_c A(X) \nodalValue{w} = \gamma \parentheses[\big]{b M(X) \nodalValue{1} - M(X; \nodalValue{u},\nodalValue{u}) \nodalValue{w}} + F_{(2)}. \end{gather*}$

We note that instead of the $ L^2$ -projection we use the interpolation of $ g$ , $f_{(1)}$ respectively $f_{(2)}$ .

Full discretization (Elliott+Styles discretization)

We perform three steps.

Given $\nodalValue{X}^n,\, \nodalValue{u}^n,\, \nodalValue{w}^n$ solve for $\nodalValue{X}^{n+1}$
$\begin{equation*} \parentheses[\big]{M_3^n + (\alpha + \varepsilon\tau) A_3^n} \nodalValue{X}^{n+1} = (M_3^n + \alpha A_3^n) \nodalValue{X}^n + \tau \delta M_3^n(\nodalValue{u}^n, \nodalValue{\surfaceNormal}) + \tau G^n, \end{equation*}$
where $\nodalValue{\surfaceNormal}^n$ is elementwise normal.
Given $\nodalValue{X}^{n+1},\, \nodalValue{u}^n,\, \nodalValue{w}^n$ solve for $\nodalValue{u}^{n+1}$
$\begin{equation*} (1+ \tau \gamma) (M\nodalValue{u})^{n+1} + \tau (A \nodalValue{u})^{n+1} - \tau \gamma M^{n+1}(\nodalValue{u}^n, \nodalValue{w}^n) \nodalValue{u}^{n+1} = (M\nodalValue{u})^n + \tau \gamma a M^{n+1} \nodalValue{1} + \tau F^{n+1}_{(1)}, \end{equation*}$
where a tensor is, namely $M_{ijkl} = \int \chi_i \chi_j \chi_k \chi_l$ and $\nodalValue{1}$ means the constant finite element function.
Given $\nodalValue{X}^{n+1},\, \nodalValue{u}^{n+1},\, \nodalValue{w}^n$ solve for $\nodalValue{w}^{n+1}$
$\begin{equation*} (M\nodalValue{w})^{n+1} + \tau D_c (A \nodalValue{w})^{n+1} + \tau \gamma M^{n+1}(\nodalValue{u}^{n+1}, \nodalValue{u}^{n+1}) \nodalValue{w}^{n+1} = (M \nodalValue{w})^n + \tau \gamma b M^{n+1} \nodalValue{1} + \tau F^{n+1}_{(2)}. \end{equation*}$

Definition at line 147 of file brusselator_algo.h.

Constructor & Destructor Documentation

Brusselator_scheme::Brusselator_scheme	(	int	argc,
		char **	argv,
		const std::string &	parameter_fname
	)

The constructor that also performs the first part before the loop enters.

Parameters

	argc	`argc` from `main`
[in]	argv	`argv` from `main`
	parameter_fname	Preferable absolute path to parameter file.

Warning: Absolute path differs on different operating systems.

Todo:: Test if folder for tmp file exists.

Definition at line 71 of file brusselator_algo.cpp.

Member Function Documentation

template<class It1 , class It2 >

void Esfem::Brusselator_scheme::calculate_velocity	(	It1	xn_first,
		It1	xn_last,
		It1	xo_first,
		It2	v_first
	)

inlineprivate

Calculate the velocity via simple differential quotient.

Parameters

[in]	xn_first	Iterator to first point of the new surface
[in]	xn_last	Iterator to last point of the new surface
[in]	xo_first	Iterator to first point of the previous surface
[out]	v_first	Iterator to the first value of the velocity

Precondition: Both xo_first and v_first point to as many elements as xn_first does.

Definition at line 438 of file brusselator_algo.h.

void Brusselator_scheme::eoc_logisticSphere ( )

Experiment, where right-hand side is calculated for a known solution.

As exact solution for the scalar valued surface equation I chose

$\begin{equation*} u(x,y,z,t) = x y e^{-6t} \quad \text{and}\quad w(x,y,z,t) = y z e^{-6t}. \end{equation*}$

For the surface evolution I chose

$\begin{equation*} \Phi(x,t) := r(t) x, \quad r(t) := \frac{r_{end} r_0}{r_{end} e^{-kt} + r_0 (1-e^{-kt})}, \end{equation*}$

where $ r(t)$ is the logistic growth function, which satisfies the following ODE

$\begin{equation*} \dot{r} = k \Bigl( 1 - \frac{r}{r_{end}}\Bigr) r,\quad r(0) = r_0. \end{equation*}$

From this it follows easily that the velocity is given by

$\begin{equation*} v(x,t) = k \Bigl(1 - \frac{r}{r_{end}}\Bigr) x. \end{equation*}$

Definition at line 116 of file brusselator_algo.cpp.

void Brusselator_scheme::eoc_mcf ( )

Dziuk mean curvature flow ESFEM.

(This is not any more true) Experiment reads as follows: Stationary surface $\surface = S^2$ with exact solution $X=e^{-t}(xy, yz, xz)$ , with the formula

$\Delta f = \Delta_{\R^{n+1}} f - H D_{\R^{n+1}} f \cdot n - D^2_{\R^{n+1}} f(n,n),$

and

$H = div(n) = \frac{n}{|x|}$

one easily sees

$\Delta(xy) = -\frac{2n+2}{|x|^2} xy,$

which implies that the right-hand side of the surface PDE vanishes.

Precondition: delta == 0, epsilon == 1, alpha == 0.

Definition at line 169 of file brusselator_algo.cpp.

void Brusselator_scheme::eoc_sdp ( )

EOC experiment for solution driven paper.

As exact solution for the surface evolution I choose

$\begin{equation*} \Phi(x,t) := r(t) x, \quad r(t) := \frac{r_{end} r_0}{r_{end} e^{-kt} + r_0 (1-e^{-kt})}. \end{equation*}$

$r(t)$ satiesfies the ODE

$\dot{r}(t) = k \left( 1 - \frac{r(t)}{r_{end}}\right) r(t),$

which implies for the velocity

$v(x,t) = k \left( 1 - \frac{r(t)}{r_{end}} \right)x =: \tilde{a} x$

As exact solution for the scalar diffusion equation I choose

$u(x,t) := x_1 x_2 e^{-6t}.$

The coupled PDE equations reads as

$\begin{gather*} \partial^{\bullet} u + u \nabla\cdot v - \Delta u = f \\ v - \Delta v - \Delta X - \delta u = g, \end{gather*}$

where $g$ is given via

$g = \left( \tilde{a} + \frac{(\alpha \tilde{a} + \varepsilon)H - \delta u}{|x|}\right) x,$

and the complicated $f$ is computed via Sage.

Precondition: Parameter $\gamma$ has to be zero.

See also: eoc_sls(), Esfem::SecOrd_op::vRhs::new_sls(), Esfem::SecOrd_op::sRhs::new_sdp_u(), Esfem::SecOrd_op::vIdata::new_sls(), Esfem::SecOrd_op::sIdata::new_1ssef()

Definition at line 242 of file brusselator_algo.cpp.

void Brusselator_scheme::eoc_sls ( )

Surface logistic sphere.

As exact solution for the surface evolution I choose

$\begin{equation*} \Phi(x,t) := r(t) x, \quad r(t) := \frac{r_{end} r_0}{r_{end} e^{-kt} + r_0 (1-e^{-kt})}. \end{equation*}$

$r(t)$ satiesfies the ODE

$\dot{r}(t) = k \left( 1 - \frac{r(t)}{r_{end}}\right) r(t),$

which implies for the velocity

$v(x,t) = k \left( 1 - \frac{r(t)}{r_{end}} \right)x =: \tilde{a} x$

I do not consider coupling. For the surface PDE I choose

$v - \alpha \Delta v - \varepsilon \Delta x - \delta u \vec{n} = f,$

where $f$ must be

$f = \left( \tilde{a} + \frac{(\alpha \tilde{a} + \varepsilon)H - \delta u}{|x|}\right) x,$

where we have used $\Delta x = - H \vec{n}$ , where $H$ is the mean curvature (without aritmetic mean) and $\vec{n}$ is the outwards pointing normal, and $\vec{n} = \frac{x}{|x|}$ . It holds $H = \frac{n}{|x|}$ , where $n$ is the dimension of the sphere. Note that $|x| = r(t)$ on the exact surface.