Commit d5b0265d by Jigyasa Watwani

1D/2D/3D solutions for all three growth models. Problem in the linear solution.

parents 28c84ab6 8b2e5a8c
......@@ -47,34 +47,37 @@ class GrowthDiffusion(object):
+ self.Dc * df.inner(df.nabla_grad(self.c), df.nabla_grad(tc))
- df.inner(self.reaction_rate * self.c, tc) ) * df.dx
def solve(self):
times = np.arange(0, self.maxtime+self.timestep, self.timestep)
fname = params['timestamp'] + '_concentration'
cFile = df.XDMFFile(fname+ '.xdmf')
savesteps = int(self.savetime/self.timestep)
maxsteps = int(self.maxtime/self.timestep)
self.time = 0.0
cFile = df.XDMFFile('%s_concentration.xdmf' %params['timestamp'])
# initial condition
r2 = "+".join(['x['+str(i)+']*x['+str(i)+']' for i in range(self.dimension)])
c0 = '1 + 0.1*cos(pi*m*(%s)/L)' % r2
c0 = df.Expression(c0, pi=np.pi, m=1, L=self.system_size, degree=1)
self.c0.assign(df.project(c0, self.SFS))
# save data
cFile.write_checkpoint(self.c0, fname, times[0])
cFile.write_checkpoint(self.c0, 'concentration', self.time)
# time stepping
for i in progressbar.progressbar(range(1, len(times))):
for steps in progressbar.progressbar(range(1, maxsteps + 1)):
# get velocity
self.sigma.t = times[i-1]
self.sigma.t = self.time
self.velocity.assign(df.project(self.sigma*self.growth_direction, self.VFS))
# solve
df.solve(self.form == 0, self.c)
# update
self.c0.assign(self.c)
# save data
cFile.write_checkpoint(self.c0, fname, times[i], append=True)
if steps % savesteps == 0:
cFile.write_checkpoint(self.c0, 'concentration', self.time, append=True)
# move mesh
displacement = df.project(self.velocity*self.timestep, self.VFS)
df.ALE.move(self.mesh, displacement)
self.time += self.timestep
cFile.close()
......@@ -95,5 +98,5 @@ if __name__ == '__main__':
gd = GrowthDiffusion(params)
gd.solve()
with open(params['timestamp'] + '_parmeters.json', "w") as fp:
with open(params['timestamp'] + '_parameters.json', "w") as fp:
json.dump(params, fp, indent=4)
\ No newline at end of file
\documentclass[11pt, reqno]{report}
\usepackage{amssymb, amsmath, datetime, xcolor, soul, graphicx}
\documentclass[11pt, reqno]{amsart}
\usepackage{amssymb,amsmath,datetime,soul,graphicx,geometry}
\usepackage[dvipsnames]{xcolor}
\linespread{1.25}
\usepackage[colorlinks=true, citecolor=blue, urlcolor=blue]{hyperref}
\usepackage[colorlinks=true,citecolor=RoyalBlue,urlcolor=RoyalBlue]{hyperref}
% shortcuts
\newcommand{\vect}[1]{\boldsymbol{#1}}
......@@ -13,243 +14,201 @@
\begin{document}
\title{Diffusion on a growing domains}
\title{Diffusion in isotropically growing domains}
\date{\currenttime, \today}
\maketitle
Largely following \href{https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0117949}{Exact Solutions of Linear Reaction-Diffusion Processes on a Uniformly Growing Domain: Criteria for Successful Colonization}, M J Simpson, PLoS One (2015).
\section{Diffusion on a growing disk}
\subsection{Kinematics}
We consider a linear reaction-diffusion process on a growing disk labelled by points $\vect{r} = (r, \varphi)$ with $0 < r < R(t)$ and $0 \leq \varphi < 2\pi$ where $R(t)$ is the increasing radius of the domain. Domain growth is associated with a \emph{radial} velocity field $\vect{v} = v(r,t) \, \uvec{r}$ which causes points on a circle of radius $r$ to move to a circle of radius $r + v(r,t) \, \tau $ in a short time $\tau$. Thus,
\eqn{r \to r^{\prime} = r + v(r, t)\tau.\label{eq:rprime}}
And, similarly,
\eqn{r + \Delta r \to r^{\prime} + \Delta r^{\prime} = r + \Delta r + v(r + \Delta r, t)\tau .\label{eq:rprimeplusdeltarprime}}
Subtracting Eq(\ref{eq:rprime}) from Eq(\ref{eq:rprimeplusdeltarprime}),
\begin{equation}
\begin{split}
\Delta r^{\prime} &= \Delta r + \tau[v(r + \Delta r, t) - v(r,t)] \\
&\approx \Delta r \left( 1 + \tau \frac{\partial v}{\partial r}\right)
\end{split}
\end{equation}
Integrating both sides,
\eq{
R^{\prime}(t) = R(t) + \tau \int_0^{R(t)} \frac{\partial v}{\partial r} dr
\section{Kinematics}
We consider a linear reaction-diffusion process on a growing $n-$dimensional domain labelled by points $\vect{r} = (r, \uvec{n})$ where $0 < r < R(t)$ with $R(t)$ is the increasing radius of the domain and $\uvec{n}$ is a unit-vector indicating a direction ($\uvec{n}$ only exists for $n \geq 2$). Domain growth is associated with a \emph{radial} velocity field $\vect{v} = v(r,t) \, \uvec{r}$ which causes points at a distance $r$ from the origin to move to a distance $r + v(r,t) \, \Delta t $ in a short time $\Delta t$. Thus points at distances $r$ and $r + \Delta r$, respectively, move to
\eqn{r \to r^{\prime} &= r + v(r, t) \;\tau,
\\
r + \Delta r \to r^{\prime} + \Delta r^{\prime} &= r + \Delta r + v(r + \Delta r, t) \; \Delta t.}
Thus
\eqn{
\Delta r^{\prime} = \Delta r + \Delta t \; \big[ v(r + \Delta r, t) - v(r,t) \big]
\approx
\Delta r \left( 1 + \Delta t \frac{\partial v}{\partial r}\right)
}
or
implying
\eqn{
R^{\prime}(t) = R(t) + \Delta t \int_0^{R(t)} \frac{\partial v}{\partial r} dr
\qquad
\implies
\qquad
\frac{dR(t)}{dt} = \int_0^{R(t)} \frac{\partial v}{\partial r} dr
.
\label{eq:growing_radius}
}
We consider uniform growth, i.e., $\partial_r v$ is independent of $r$ but potentially depends on time $t$, so that we have $\partial_r v = \sigma(t)$. Combining this with \eqref{eq:growing_radius} gives
We consider uniform growth, i.e., $\partial_r v \equiv \sigma(t)$ is independent of $r$ but could potentially depend on time $t$. Thus
\eqn{
\frac{\partial v}{\partial r} = \sigma(t) = \frac{1}{R(t)} \frac{dR(t)}{dt}.
\sigma(t) = \frac{1}{R(t)} \frac{dR(t)}{dt}.
\label{eq:growth_rate}
}
Assuming that the circular domain elongates with its center fixed, i.e., $v(0,t) = 0$, integrating \eqref{eq:growth_rate} gives
We impose the boundary condition $v(0,t) = 0$ using which, on integrating \eqref{eq:growth_rate}, we find
\eqn{
v(r,t) = \frac{r}{R(t)} \frac{dR(t)}{dt}.
v(r,t) = \frac{r}{R(t)} \frac{dR(t)}{dt} = \sigma(t) \; r.
}
\subsection{Mass conservation}
\section{Mass conservation}
We now consider the conservation of a mass density $C(\vect{r}, t)$, assuming that it evolves according to a linear reaction-diffusion mechanism. The associated conservation statement on the growing domain can be written as
\eqn{
\partial_t C = D \nabla^2 C - \nabla \cdot (\vect{v} C) + k \, C,
\label{eq:general_reaction_diffusion}
}
where $D>0$ and $k$ are the diffusion constant and reaction rate respectively.
We assume azimuthal symmetry $C(\vect{r}, t) = C(r, t)$, i.e., the dynamics is independent of $\varphi$. Then
where $D>0$ and $k$ are the diffusion constant and reaction rate respectively. Assuming a solution that only depends on the growth direction, i.e., $C(\vect{r}, t) = C(r, t)$, we get
\eqn{
\frac{\partial C}{\partial t} = D \, \frac{1}{r} \, \frac{\partial}{\partial r} \left(r \frac{\partial C}{\partial r} \right) - \frac{1}{r} \, \frac{\partial (r v C)}{\partial r} + k\,C,
\frac{\partial C}{\partial t} = \frac{D}{r^{n-1}} \, \frac{\partial}{\partial r} \left(r^{n-1} \frac{\partial C}{\partial r} \right) - \frac{1}{r^{n-1}} \, \frac{\partial (r^{n-1} v C)}{\partial r} + k\,C,
\\
= \frac{D}{r^{n-1}} \, \frac{\partial}{\partial r} \left(r^{n-1} \frac{\partial C}{\partial r} \right) - \frac{\sigma(t)}{r^{n-1}} \, \frac{\partial (r^{n} C)}{\partial r} + k\,C.
\label{eq:azimuthal_symmetry}
}
and impose zero diffusive flux conditions $\partial_r C = 0$ at $r=0$ and at $r=R(t)$. With some given initial conditions $C(r, 0)$ our aim is to find the exact solution of \eqref{eq:azimuthal_symmetry}.
We impose zero diffusive flux conditions $\partial_r C = 0$ at $r=0$ and at $r=R(t)$. With some given initial conditions $C(r, 0)$ our aim is to find the exact solution of \eqref{eq:azimuthal_symmetry}.
\subsection{Exact solution}
\section{Rescaling}
We first transform the spatial variable to a fixed domain $\xi = \dfrac{r}{R(t)}$. Then
We first transform the spatial variable to a fixed domain:
\eqn{ r \to \xi = \frac{r}{R}.}
This leads to
\eqn{
\frac{\partial C}{\partial t} \to \frac{\partial C}{\partial t} - \frac{v}{R(t)} \frac{\partial C}{\partial \xi}
\frac{\partial}{\partial t} \to \frac{\partial}{\partial t} - \sigma \; \xi \frac{\partial}{\partial \xi},
\qquad \qquad
\frac{\partial}{\partial r} \to \frac{1}{R} \frac{\partial}{\partial \xi}
}
which implies
\eqn{
\frac{\partial C}{\partial r} \to \frac{1}{R(t)} \frac{\partial C}{\partial \xi}
\frac{\partial C}{\partial t} =
\frac{D }{\xi^{n-1} R^2} \, \frac{\partial}{\partial \xi} \left(\xi^{n-1} \frac{\partial C}{\partial \xi} \right) + (k - n \, \sigma) C.
}
Then Eq(\ref{eq:azimuthal_symmetry}) becomes:
\eqn{\frac{\partial C}{\partial t} = \frac{D}{\xi R^2(t)} \frac{\partial}{\partial \xi} \left( \xi \frac{\partial C}{\partial \xi} \right) + \left(k - 2 \sigma(t)\right) C}
We now transform the variable $t: t \to T = \int_0^t \frac{D}{R^2(s)} ds$, giving
\eqn{\frac{\partial C}{\partial T} = \frac{1}{\xi} \frac{\partial}{\partial \xi} \left( \xi \frac{\partial C}{\partial \xi} \right) + \frac{R^2(t)}{D}\left(k - 2 \sigma(t) \right) C}
Finally, we express the second term on the RHS as a function of $T$ to write:
\eqn{\frac{\partial C}{\partial T} = \frac{1}{\xi} \frac{\partial}{\partial \xi} \left( \xi \frac{\partial C}{\partial \xi} \right) + f(T) C}
The above equation can be solved easily by separation of variables. Writing $C(\xi, T) = u(\xi) w(T)$,
\eqn{\frac{1}{w} \frac{dw}{dt} = \frac{1}{u \xi} \frac{d}{d \xi} \left( \xi \frac{du}{d \xi} \right) + f(T)}
Thus,
\eqn{\frac{1}{w} \frac{dw}{dT} - f(T) = \frac{1}{u \xi} \frac{d}{d \xi} \left( \xi \frac{du}{d \xi} \right) = -n^2}
Solving the $w$ equation first:
\eqn{w(T) = w(0) \exp\left[-n^2 T + \int_0^T f(T^{\star}) dT^\star \right]}
The $u$ equation reads:
\eqn{\xi \frac{d^2 u}{d \xi^2} + \frac{du}{d \xi} + n^2 u \xi = 0}
This is the Bessel equation of zeroth order, having solutions
\eqn{u(\xi) = c_1 J_0 (\xi n) + c_2 yY_0 (\xi n)}
Thus, the exact solution in the $(\xi, T)$ coordinates is :
\eqn{C(\xi, T) = \left[ a_1 J_0 (\xi n) + a_2 Y_0 (\xi n)\right] \exp\left[-n^2 T + \int_0^T f(T^\star) dT^\star \right]}
Since we want the solution to be bounded at $\xi = 0$, we keep only the $J_0$ term:
\eqn{C(\xi, T) = a_1 J_0 (\xi n) \exp\left[-n^2 T + \int_0^T f(T^\star) dT^\star \right]}
If we had Dirichlet boundary conditions at $\xi = 1$, $n$ would be the solution of
\eqn{J_0 (n) = 0}
i.e $n$'s are the roots of the zeroth bessel function.
Thus, the most general solution is:
\eqn{C(\xi, T) = \sum_{n} a_n J_0(\eta_n \xi) \exp\left[-\eta_n^2 T + \int_0^T f(T^\star) dT^\star \right]
\label{general_sol}}
where $\eta_n$ is the $n$th root of the zeroth Bessel function for Dirichlet boundary conditions. \\$a_n$ is determined by the initial conditions. Let us take
\eqn{C(\xi, 0) = J_0(\eta_m \xi) = \sum_{n} a_n J_0(\eta_n \xi)}
where the last equality follows by putting $T=0$ in Eq(\ref{general_sol}). Now using the orthogonality relation
\eq{\int_0^a J_v \left( \eta_{\nu p}\frac{\xi}{a}\right)J_v \left( \eta_{\nu q}\frac{\xi}{a}\right) \xi d \xi = \frac{a^2}{2}[J_{\nu+1} (\eta_{\nu p})]^2 \delta_{pq}}
where $\eta_{\nu p}$ and $\eta_{\nu q}$ are the $p$th and $q$th roots of the $\nu$th bessel function, we get
\eqn{a_{n=m} = 1, a_{n \neq m} = 0}.
In summary, the solution to our problem for the initial condition $C(\xi, 0) = J_0(\eta_m \xi)$ where $\eta_m$ is the mth root of $j_0$ is:
\eqn{C(\xi, T) = J_0(\eta_{m} \xi) \exp \left[-\eta_m^2 T + \int_0^T f(T^\star)dT^\star \right]}.
\subsection{Specific examples}
\subsubsection*{No growth}
In this case, $R(t) = R, \sigma(t) = 0, f(T^\star) = R^2 k/D$.
This gives \eq{C(\xi, T) = J_0 (\eta_m \xi) e^{-\eta_m^2 T} e^{R^2 k T/D}}
Plugging in $T = D t/R^2 $ and $\xi = r/R$,
\eqn{C(r,t) = J_0 (\eta_m r/R) e^{-\frac{\eta_m^2 D t}{R^2}}e^{kt}}
\subsubsection{Exponential domain growth}
In this case, $R(t) = R_0 e^{\alpha t}, \sigma(t)=\alpha$
\eq{ T = \frac{D(1 - e^{-2 \alpha t})}{2 \alpha R_0^2}}
And \eq{f( T^\star) = \frac{R_0^2 (k - 2 \alpha)}{D - 2 \alpha R_0^2 T^{\star}}}
Then,
\eq{C(\xi, T) = J_0(\eta_m \xi) e^{-\eta_m^2 T} \left(\frac{D}{D - 2\alpha R_0^2 T}\right)^{(k-2\alpha)/2\alpha}}
and \eqn{C(r,t) = J_0\left(\eta_m \frac{r}{R(t)}\right) \exp\left( -\frac{\eta_m^2 D(1 - e^{-2 \alpha t})}{2 \alpha R_0^2}\right) e^{(k - 2\alpha)t}}
\subsubsection{Linear domain growth}
In this case, $R(t)=R_0 + bt, \sigma(t) = b/R(t)$.\\
Further, \eq{T = \frac{D t}{R(t) R_0}}, or,
\eq{t = \frac{R_0^2 T}{D - R_0 b T}}
And \eq{f(T^\star) = \frac{R_0}{(D - R_0 bT^\star)^2} (k D R_0 - 2 b D + 2bR_0T^\star)}
Then
\eq{C(\xi, T) = J_0(\eta_m \xi) e^{-\eta_m^2 T} e^{\frac{(kR_0^2 - 2b + 2 R_0)T}{D - R_0 bT}}\left(\frac{D - R_0 bT}{D}\right)^\frac{2 R_0}{b}}
\eqn{C(r, t) = J_0\left(\eta_{m} \frac{n r}{R(t)}\right)e^{-\frac{\eta_m^2 D t}{R(t) R_0}} \left(\frac{R_0}{R(t)}\right)^{2R_0/b} e^{kt} e^{\frac{2( R_0-b)t}{R_0^2}}}
\section{Diffusion in a growing sphere}
\subsection{Kinematics}
We consider a linear reaction-diffusion process on a growing sphere labelled by points $\vect{r} = (r, \theta, \varphi)$ with $0 < r < R(t)$, $0 \leq \varphi < 2\pi$ and $0 \leq \theta \leq \pi$ where $R(t)$ is the increasing radius of the domain. Domain growth is associated with a \emph{radial} velocity field $\vect{v} = v(r,t) \, \uvec{r}$ which causes points on the surface of the sphere of radius $r$ to move to a circle of radius $r + v(r,t) \, \tau $ in a short time $\tau$. Thus,
\eqn{r \to r^{\prime} = r + v(r, t)\tau.\label{eq:rprime}}
And, similarly,
\eqn{r + \Delta r \to r^{\prime} + \Delta r^{\prime} = r + \Delta r + v(r + \Delta r, t)\tau .\label{eq:rprimeplusdeltarprime}}
Subtracting Eq(\ref{eq:rprime}) from Eq(\ref{eq:rprimeplusdeltarprime}),
\begin{equation}
\begin{split}
\Delta r^{\prime} &= \Delta r + \tau[v(r + \Delta r, t) - v(r,t)] \\
&\approx \Delta r \left( 1 + \tau \frac{\partial v}{\partial r}\right)
\end{split}
\end{equation}
Integrating both sides,
\eq{
R^{\prime}(t) = R(t) + \tau \int_0^{R(t)} \frac{\partial v}{\partial r} dr
}
or
Next, we transform the variable:
\eqn{
\frac{dR(t)}{dt} = \int_0^{R(t)} \frac{\partial v}{\partial r} dr
.
\label{eq:growing_radius}
t \to \tau = \int_0^t ds \; \frac{D}{R^2(s)}.
}
We consider uniform growth, i.e., $\partial_r v$ is independent of $r$ but potentially depends on time $t$, so that we have $\partial_r v = \sigma(t)$. Combining this with \eqref{eq:growing_radius} gives
Then
\eqn{
\frac{\partial v}{\partial r} = \sigma(t) = \frac{1}{R(t)} \frac{dR(t)}{dt}.
\label{eq:growth_rate}
\frac{\partial C}{\partial \tau} =
\frac{1}{\xi^{n-1}} \, \frac{\partial}{\partial \xi} \left(\xi^{n-1} \frac{\partial C}{\partial \xi} \right) + f(\tau) C
\label{eq:nondim_eqn}
}
Assuming that the circular domain elongates with its center fixed, i.e., $v(0,t) = 0$, integrating \eqref{eq:growth_rate} gives
where
\eqn{
v(r,t) = \frac{r}{R(t)} \frac{dR(t)}{dt}.
f(t(\tau)) = \frac{R(t(\tau))^2}{D} \big[k - n \, \sigma(t(\tau)) \big].
}
\subsection{Mass conservation}
\section{Exact solutions}
We now consider the conservation of a mass density $C(\vect{r}, t)$, assuming that it evolves according to a linear reaction-diffusion mechanism. The associated conservation statement on the growing domain can be written as
We use separation of variables to solve \eqref{eq:nondim_eqn}. Writing $C(\xi, \tau) = u(\xi) \; w(\tau)$,
\eqn{\frac{1}{w} \frac{dw}{d\tau} - f(\tau)= \frac{1}{u \, \xi^{n-1}} \frac{d}{d \xi} \left( \xi^{n-1} \frac{du}{d \xi} \right).}
With $-\omega^2$ as a separation constant, we get
\eqn{
\partial_t C = D \nabla^2 C - \nabla \cdot (\vect{v} C) + k \, C,
\label{eq:general_reaction_diffusion}
\frac{dw}{d\tau} = \big[ f(\tau) - \omega^2 \big] w
\qquad
\implies
\qquad
w(\tau) = w(0) \exp\left[-\omega^2 \tau + \int_0^\tau ds\; f(s) \right],
}
where $D>0$ and $k$ are the diffusion constant and reaction rate respectively.
We assume azimuthal symmetry $C(\vect{r}, t) = C(r, t)$, i.e., the dynamics is independent of $\varphi$. Then
and
\eqn{
\frac{\partial C}{\partial t} = D \, \frac{1}{r^2} \, \frac{\partial}{\partial r} \left(r^2 \frac{\partial C}{\partial r} \right) - \frac{1}{r^2} \, \frac{\partial (r^2 v C)}{\partial r} + k \, C,
\label{eq:azimuthal_symmetry}
}
and impose zero diffusive flux conditions $\partial_r C = 0$ at $r=0$ and at $r=R(t)$. With some given initial conditions $C(r, 0)$ our aim is to find the exact solution of \eqref{eq:azimuthal_symmetry}.
\subsection{Exact solution}
We first transform the spatial variable to a fixed domain $\xi = \dfrac{r}{R(t)}$. Then
\frac{d^2 u}{d\xi^2} + \frac{(n-1)}{\xi} \frac{du}{d\xi} + \omega^2 u = 0.
\label{eq:u_eqn}
}
We need to find solutions of \eqref{eq:u_eqn} with the boundary conditions that $u'(\xi) = 0$ at both $\xi=0$ and $\xi=1$.
\begin{itemize}
\item When $n=1$, we get
\eqn{u'' + \omega^2 u = 0,}
and thus we need $\omega=m\pi, m \in \mathbb{Z}$ to satisfy the boundary conditions. The general solution is thus
\eqn{u(\xi, \tau) = \sum_{m=0}^{\infty} a_m \cos\left(m \pi \xi \right) \; \exp\left[-m^2 \pi^2 \tau + \int_0^\tau ds\; f(s) \right].}
%
\item When $n=2$, we get
\eqn{\xi \, u'' + u' + \omega^2 \, \xi \, u = 0,}
the solutions of which are the zeroth-order Bessel functions $J_0(\omega\xi)$ and $Y_0(\omega\xi)$. Requiring $u$ to be bounded at $\xi=0$, we neglect the $Y_0$ solution. Imposing the boundary condition
\eqn{\frac{dJ_0(z)}{dz}\bigg\vert_{z=\omega}=0 \qquad \implies J_1(\omega)=0,
\label{eq:bc_n2}}
determines the allowed values of $\omega$ as the zeroes of the 1st-order Bessel function. The general solution is thus
\eqn{u(\xi, \tau) = \sum_{\omega} a_{\omega} \, J_0(\omega \xi) \; \exp\left[-\omega^2 \tau + \int_0^\tau ds\; f(s) \right],}
where the summation is over those values of $\omega$ satisfying \eqref{eq:bc_n2}.
%
\item When $n=3$, we get
\eqn{\xi \, u'' + 2 \, u' + \omega^2 \, \xi \, u = 0,}
the solutions of which are the zeroth-order \emph{spherical} Bessel functions
\eqn{
\frac{\partial C}{\partial t} \to \frac{\partial C}{\partial t} - \frac{v}{R(t)} \frac{\partial C}{\partial \xi}
}
j_0(\omega \xi) = \frac{\sin (\omega \xi)}{\omega \xi},
\qquad \mathrm{and} \qquad
y_0(\omega \xi) = -\frac{\cos (\omega \xi)}{\omega \xi}.}
Requiring $u$ to be bounded at $\xi=0$, we neglect the $y_0$ solution. Imposing the boundary condition
\eqn{\frac{dj_0(z)}{dz}\bigg\vert_{z=\omega}=0 \qquad \implies \qquad
\tan \omega = \omega,
\label{eq:bc_n3}}
determines the allowed values of $\omega$. The general solution is thus
\eqn{u(\xi, \tau) = \sum_{\omega} a_{\omega} \, \frac{\sin (\omega \xi)}{\omega \xi} \; \exp\left[-\omega^2 \tau + \int_0^\tau ds\; f(s) \right],}
where the summation is over those values of $\omega$ satisfying \eqref{eq:bc_n3}.
\end{itemize}
\section{Three cases}
\subsection{Non-growing domain}
With $R(t) = R_0$ in this case, we get
\eqn{\sigma(t)=0, \qquad \tau=\frac{Dt}{R_0^2}, \qquad f=\frac{R_0^2k}{D}.}
Hence $\tau \to\infty$ as $t\to\infty$. Thus
\begin{itemize}
%
\item $n=1$
\eqn{
\frac{\partial C}{\partial r} \to \frac{1}{R(t)} \frac{\partial C}{\partial \xi}
}
Then Eq(\ref{eq:azimuthal_symmetry}) becomes:
\eqn{\frac{\partial C}{\partial t} = \frac{D}{\xi^2 R^2(t)} \frac{\partial}{\partial \xi} \left( \xi^2 \frac{\partial C}{\partial \xi} \right) + \left(k - 3 \sigma(t)\right) C}
We now transform the variable $t: t \to T = \int_0^t \frac{D}{R^2(s)} ds$, giving
\eqn{\frac{\partial C}{\partial T} = \frac{1}{\xi^2} \frac{\partial}{\partial \xi} \left( \xi^2 \frac{\partial C}{\partial \xi} \right) + \frac{R^2(t)}{D}\left(k - 3 \sigma(t) \right) C}
Finally, we express the second term on the RHS as a function of $T$ to write:
\eqn{\frac{\partial C}{\partial T} = \frac{1}{\xi^2} \frac{\partial}{\partial \xi} \left( \xi^2 \frac{\partial C}{\partial \xi} \right) + f(T) C}
The above equation can be solved easily by separation of variables. Writing $C(\xi, T) = u(\xi) v(T)$,
\eqn{\frac{1}{v} \frac{dv}{dt} = \frac{1}{u \xi^2} \frac{d}{d \xi} \left( \xi^2 \frac{du}{d \xi} \right) + f(T)}
Thus,
\eqn{\frac{1}{v} \frac{dv}{dT} - f(T) = \frac{1}{u \xi^2} \frac{d}{d \xi} \left( \xi^2 \frac{du}{d \xi} \right) = -n^2}
Solving the $v$ equation first:
\eqn{v(T) = v(0) \exp\left[-n^2 T + \int_0^T f(T^{\star}) dT^\star \right]}
And the $u$ equation reads:
\eqn{\frac{d^2 u}{d \xi^2} + \frac{2}{\xi} \frac{du}{d \xi} + n^2 u = 0}
This is the Bessel equation of zeroth order, having solutions
\eqn{u(\xi) = c_1 j_0 (n \xi) + c_2 y_0 (n \xi)}
where $j_0, y_0$ are spherical bessel functions given by
\eq{j_0 (n \xi) = \frac{\sin(n \xi)}{n \xi}}
\eq{y_0 (n \xi) = -\frac{\cos(n \xi)}{n \xi}}
Thus, the exact solution in the $(\xi, T)$ coordinates is :
\eqn{C(\xi, T) = \left[ a_1 j_0 (n \xi) + a_2 y_0 (n \xi)\right] \exp\left[-n^2 T + \int_0^T f(T^\star) dT^\star \right]}
Since we want the solution to be bounded at $\xi = 0$, we keep only the $j_0$ term:
\eqn{C(\xi, T) = a_1 j_0 (n \xi) \exp\left[-n^2 T + \int_0^T f(T^\star) dT^\star \right]}
If we have Dirichlet boundary conditions at $\xi =0, 1$ then $n$ can be determined to be the solution of $j_0(n) = 0$, i.e.
\eq{n = m \pi, m \in \mathbb{Z}}
Thus, the most general solution is:
\eq{C(\xi, T) = \sum_{m} a_m j_0(m \pi \xi) \exp\left[-(m \pi)^2 T + \int_0^T f(T^\star) dT^\star \right]}
We now determine $a_m$ by the initial conditions. Let us take
\eqn{C(\xi, 0) = j_0(p \pi \xi) = \sum_m a_m j_0(m \pi \xi)}
Using the orthogonality relation
\eqn{\int_0^a j_{\nu}\left(\eta_{\nu p}\frac{\xi}{a}\right) j_{\nu}\left(\eta_{\nu q}\frac{\xi}{a}\right) \xi^2 d \xi = \frac{a^3}{2} [j_{\nu + 1}(\eta_{\nu p})]^2 \delta_{pq}}
where $\eta_{\nu p}$ and $\eta_{\nu q}$ are the $p$th and $q$th roots of the $\nu$th bessel function, we get
\eqn{a_{m=p}=1, a_{m \neq p} = 0}
And thus the solution for initial condition $C(\xi, 0) = j_0(p \pi \xi)$ and Dirichlet boundaries at $\xi = 0,1$ is:
\eqn{C(\xi, T) = j_0(p \pi \xi) \exp\left[-(p \pi)^2 T + \int_0^T f(T^\star) dT^\star\right]}
\subsection{Specific examples}
\subsubsection*{No growth}
In this case, $R(t) = R, \sigma(t) = 0, f(T^\star) = R^2 k/D$.
This gives \eq{C(\xi, T) = j_0(p \pi \xi) e^{-p^2 \pi^2 T} e^{R^2 k T/D}}
Plugging in $T = D t/R^2 $ and $\xi = r/R$,
\eqn{C(r,t) = j_0\left(p \pi \frac{r}{R(t)}\right) e^{-\frac{p^2 \pi^2 D t}{R^2}}e^{kt}}
\subsubsection*{Exponential domain growth}
In this case, $R(t) = R_0 e^{\alpha t}, \sigma(t)=\alpha$
\eq{ T = \frac{D(1 - e^{-2 \alpha t})}{2 \alpha R_0^2}}
And \eq{f( T^\star) = \frac{R_0^2 (k - 3 \alpha)}{D - 2 \alpha R_0^2 T^{\star}}}
Then,
\eq{C(\xi, T) = j_0(p \pi \xi) e^{-(p \pi)^2 T} \left(\frac{D}{D - 2 \alpha R_0^2 T}\right)^{(k-3\alpha)/2\alpha}}
and \eqn{C(r,t) = j_0\left(p \pi \frac{r}{R(t)}\right) \exp\left[-\frac{(p \pi)^2 D(1 - e^{-2 \alpha t})}{2 \alpha R_0^2} \right] e^{(k-3\alpha)t}}
\subsubsection*{Linear domain growth}
In this case, $R(t)=R_0 + bt, \sigma(t) = b/R(t)$.\\
Further, \eq{T = \frac{D t}{R(t) R_0}}, or,
\eq{t = \frac{R_0^2 T}{D - R_0 b T}}
And \eq{f(T^\star) = \frac{R_0}{(D - R_0 bT^\star)^2} (k D R_0 - 3 b D + 3bR_0T^\star)}
Then
\eq{C(\xi, T) = j_0(p \pi \xi) e^{-p^2 \pi^2 T} e^{\frac{(k R_0 + 3 - 3b) R_0 T}{D - R_0 b T}}\left(\frac{D - R_0 b T}{D}\right)^{3/b}}
\eqn{C(r, t) = \left(\frac{R_0}{R(t)} \right)^{3/b} j_0\left(\frac{p \pi r}{R(t)}\right)e^{-\frac{p^2 \pi^2 D t}{R(t) R_0}} e^{kt} e^{\frac{3(1-b)t}{R_0}} }
\end{document}
C(r,t) = \sum_{m=0}^{\infty} a_m \cos\left(\frac{m \pi r}{R_0} \right) \; \exp\left(-\frac{m^2 \pi^2 D t}{R_0^2} + k t \right),
\qquad \mathrm{with} \quad m \in \mathbb{Z}
}
%
\item $n=2$
\eqn{C(r, t) = \sum_{\omega} a_{\omega} \, J_0\left(\frac{\omega r}{R_0} \right) \; \; \exp\left(-\frac{\omega^2 D t}{R_0^2} + k t \right),
\qquad \mathrm{with} \quad J_1(\omega)=0.
}
%
\item $n=3$
\eqn{C(r, t) = \sum_{\omega} a_{\omega} \, \frac{R_0}{\omega r} \; \sin \left(\frac{\omega r}{R_0}\right) \; \; \exp\left(-\frac{\omega^2 D t}{R_0^2} + k t \right),
\qquad \mathrm{with} \quad \tan \omega = \omega.}
\end{itemize}
\subsection{Exponential domain growth}
With $R(t) = R_0 e^{\alpha t}$ in this case, we get,
\eqn{\sigma(t)=\alpha, \qquad \tau=\frac{D \; (1-e^{-2\alpha t})}{2\alpha R_0^2} , \qquad \int_0^\tau f(s) ds = \ln\left(\frac{D}{D - 2 \alpha R_0^2 \tau}\right)^{(k-n\alpha)/2\alpha} .}
Hence $\tau \to D/(2\alpha R_0^2)$ as $t\to\infty$. Thus
\begin{itemize}
%
\item $n=1$
\eqn{C(r,t) = \sum_{m=0}^\infty a_m \cos\left( \frac{m \pi r}{R(t)}\right) \exp\left[{-\frac{m^2 \pi^2 D}{2 \alpha R_0^2} (1 - e^{-2\alpha t})+t(k-\alpha)}\right],
\qquad \mathrm{with} \quad m \in \mathbb{Z}}
%
\item $n=2$
\eqn{C(r,t) = \sum_{\omega}^\infty a_\omega J_0\left( \frac{\omega r}{R(t)}\right) \exp\left[{-\frac{\omega^2 D}{2 \alpha R_0^2} (1 - e^{-2\alpha t})+t(k-2\alpha)}\right],
\qquad \mathrm{with} \quad J_1(\omega)=0.}
%
\item $n=3$
\eqn{C(r,t) = \sum_{\omega}^\infty a_\omega \frac{ R(t)}{\omega r} \sin\left( \frac{\omega r}{R(t)}\right) \exp\left[{-\frac{\omega^2 D}{2 \alpha R_0^2} (1 - e^{-2\alpha t})+t(k-3\alpha)}\right],
\qquad \mathrm{with} \quad \tan \omega = \omega.}
\end{itemize}
\subsection{Linear domain growth}
With $R(t) = R_0 + \alpha t$ in this case, we get
\eqn{\sigma(t) = \frac{\alpha}{R_0 + \alpha t}, \qquad\tau=\frac{D \; t}{R_0 (R_0 + \alpha t)} , \qquad \int_0^\tau f(s) ds = \frac{R_0^2 k \tau}{(D - \alpha R_0 \tau)} + \ln \left( \frac{D - R_0 \alpha \tau}{D}\right)^{n}.}
Hence $\tau \to D/(\alpha R_0)$ as $t\to\infty$. Thus
\begin{itemize}
%
\item $n=1$
\eqn{C(r,t) = \sum_{m=0}^\infty a_m \cos\left( \frac{m \pi r}{R(t)}\right) \exp\left[{-\frac{m^2 \pi^2 D t}{R(t) R_0} + tk}\right] \left(\frac{R_0}{R(t)}\right),
\qquad \mathrm{with} \quad m \in \mathbb{Z}}
%
\item $n=2$
\eqn{C(r,t) = \sum_{\omega}^\infty a_\omega J_0\left( \frac{\omega r}{R(t)}\right) \exp\left[{-\frac{\omega^2 D t}{R(t) R_0} + tk}\right] \left(\frac{R_0}{R(t)}\right)^2,
\qquad \mathrm{with} \quad J_1(\omega)=0.}
%
\item $n=3$
\eqn{C(r,t) = \sum_{\omega}^\infty a_\omega \frac{ R(t)}{\omega r} \sin\left( \frac{\omega r}{R(t)}\right) \exp\left[{-\frac{\omega^2 D t}{R(t) R_0} + tk} \right] \left(\frac{R_0}{R(t)}\right)^3,
\qquad \mathrm{with} \quad \tan \omega = \omega.}
\end{itemize}
\end{document}
{
"dimension" : 2,
"resolution" : 5,
"resolution" : 15,
"system_size" : 1,
"Dc" : 0.1,
"growth" : "exponential",
"growth" : "linear",
"growth_parameter" : 0.1,
"reaction_rate" : 0.0,
"timestep" : 0.01,
"maxtime" : 1.0
"savetime": 0.1,
"maxtime" : 5.0
}
import dolfin as df
import numpy as np
import vedo as vd
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import progressbar
import os
import json
import h5py
from tempfile import TemporaryDirectory
def get_data(params, DIR=''):
savesteps = int(params['maxtime']/params['savetime'])
times = np.arange(savesteps+1) * params['savetime']
# Read mesh geometry from h5 file
var = 'concentration'
h5 = h5py.File(os.path.join(DIR, '%s_%s.h5' % (params['timestamp'], var)), "r")
# should be in the loop if remeshing
topology = np.array(h5['%s/%s_0/mesh/topology'%(var,var)])
geometry = []
for i in range(len(times)):
geometry.append(np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,i)]))
h5.close()
geometry = np.array(geometry)
# create a mesh
if params['dimension']==1:
mesh = df.IntervalMesh(params['resolution'], 0, params['system_size'])
else:
mesh = df.Mesh()
editor = df.MeshEditor()
editor.open(mesh,
type='triangle' if params['dimension']==2 else 'tetrahedron',
tdim=params['dimension'], gdim=params['dimension'])
editor.init_vertices(len(geometry[0]))
editor.init_cells(len(topology))
for i in range(len(geometry[0])):
editor.add_vertex(i, geometry[0][i])
for i in range(len(topology)):
editor.add_cell(i, topology[i])
editor.close()
# Read concentration data
concentration = np.zeros((len(times), len(geometry[0])))
cFile = os.path.join(DIR, '%s_%s.xdmf' % (params['timestamp'],var))
with df.XDMFFile(cFile) as cFile:
for steps in range(savesteps+1):
mesh.coordinates()[:] = geometry[steps]
VFS = df.FunctionSpace(mesh, 'P', 1)
c = df.Function(VFS)
cFile.read_checkpoint(c, var, steps)
concentration[steps] = c.compute_vertex_values(mesh)
return (times, geometry, topology, concentration)
def visualize(params, DIR='', offscreen=False):
n_cmap_vals = 16
scalar_cmap = 'viridis'
(times, geometry, topology, concentration) = get_data(params, DIR)
if params['dimension']==2:
geometry = np.dstack((geometry, np.zeros(geometry.shape[0:2])))
cmin, cmax = np.min(concentration), np.max(concentration)
plotter = vd.plotter.Plotter(axes=0)
poly = vd.utils.buildPolyData(geometry[0], topology)
scalar_actor = vd.mesh.Mesh(poly)
#scalar_actor.computeNormals(points=True, cells=True)
scalar_actor.pointdata['concentration'] = concentration[0]
scalar_actor.cmap(scalar_cmap, concentration[0], vmin=cmin, vmax=cmax, n=n_cmap_vals)
scalar_actor.addScalarBar(title = r'$c$',
pos=(0.8, 0.04), nlabels=2,
titleYOffset=15, titleFontSize=28, size=(100, 600))
plotter += scalar_actor
def update(idx):
scalar_actor.points(pts=geometry[idx], transformed=False)
scalar_actor.pointdata['concentration'] = concentration[idx]
def slider_update(widget, event):
value = widget.GetRepresentation().GetValue()
idx = (abs(times-value)).argmin()
update(idx)
slider = plotter.addSlider2D(slider_update, pos=[(0.1,0.94), (0.5,0.94)],
xmin=times[0], xmax=times.max(),
value=times[0], title=r"$t/\tau$")
vd.show(interactive=(not offscreen), zoom=0.8)
# make movie
if offscreen:
FPS = 10
movFile = '%s.mov' % params['timestamp']
fps = float(FPS)
command = "ffmpeg -y -r"
options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
tmp_dir = TemporaryDirectory()
get_filename = lambda x: os.path.join(tmp_dir.name, x)
for tt in range(len(times)):
idx = (abs(times-times[tt])).argmin()
update(idx)
slider.GetRepresentation().SetValue(times[tt])
fr = get_filename("%03d.png" % tt)
vd.io.screenshot(fr)
os.system(command + " " + str(fps)
+ " -i " + tmp_dir.name + os.sep
+ "%03d.png " + options + " " + movFile)
tmp_dir.cleanup()
def visualize1(params, DIR='', offscreen=False):
n_cmap_vals = 10
scalar_cmap = 'viridis'
vector_color = '#ffffff'
vector_scale = 0.1
savesteps = int(params['maxtime']/params['savetime'])
times = params['alpha'] * np.arange(savesteps+1) * params['savetime']
halftime = int(len(times)/2)
mesh = get_mesh_from_xdmf(DIR, params)
print('Reading polarity')
polarity = np.zeros((len(times), mesh.num_vertices(), 3))
polarity_functions = []
VFS = df.VectorFunctionSpace(mesh, 'P', 1)
p = df.Function(VFS)
with df.XDMFFile(os.path.join(DIR,'%s_polarity.xdmf' % params['timestamp'])) as pFile:
for steps in progressbar.progressbar(range(savesteps+1)):
pFile.read_checkpoint(p, 'polarity', steps)
polarity_functions.append(p)
vector_values = p.compute_vertex_values(mesh)
polarity[steps] = vector_values.reshape(3, int(vector_values.shape[0]/3)).T
pmag = np.linalg.norm(polarity, axis=2)
pmin, pmax = 0, np.max(pmag)
polarity = polarity / pmax
# load previously computed convergence data
data = np.loadtxt(os.path.join(DIR, '%s.dat' % params['timestamp']))
fig1, ax = plt.subplots(1)
fig1.subplots_adjust(left=0.15, bottom=0.15, right=0.9, top=0.85)
ax.plot(times, data[:,1])
ax.set_xlabel(r'$\alpha t$')
ax.set_ylabel(r'$\int_{\Gamma} (\mathbf{p}_{t+\Delta t}-\mathbf{p}_{t})^2$')
axins = inset_axes(ax, width="50%", height="40%")
axins.plot(times[halftime:], data[halftime:,1])
fig1.savefig(os.path.join(DIR, '%s_pnorm.pdf' % params['timestamp']), transparent=True)
VSH_labels = [r'$\mathbf{\Psi}_{10}$',
r'$\mathbf{\Phi}_{10}$',
r'$\mathbf{\Psi}_{11}$',
r'$\mathbf{\Phi}_{11}$',
r'$\mathbf{\Psi}_{20}$',
r'$\mathbf{\Phi}_{20}$',
r'$\mathbf{\Psi}_{21}$',
r'$\mathbf{\Phi}_{21}$',
r'$\mathbf{\Psi}_{22}$',
r'$\mathbf{\Phi}_{22}$'
]
fig2, ax = plt.subplots(1)
fig2.subplots_adjust(left=0.15, bottom=0.15, right=0.9, top=0.85)
ax.set_xlim(times[halftime:].min(), times.max())
for i in range(len(VSH_labels)):
ax.plot(times[halftime:], data[halftime:, i+2], label=VSH_labels[i], lw=1)
ax.set_xlabel(r'$\alpha t$')
ax.set_ylabel(r'$\sqrt{p^{\Psi}_{lm} (p^{\Psi}_{lm})^{\ast}}, \ \sqrt{p^{\Phi}_{lm} (p^{\Phi}_{lm})^{\ast}}$')
ax.legend(loc=0, ncol=5, fontsize='small', bbox_to_anchor=(0.1, 1))
fig2.savefig(os.path.join(DIR, '%s_modes.pdf' % params['timestamp']), transparent=True)
if not offscreen:
fig1.show()
fig2.show()
(_, _, _, _, normal, _, _) = compute_geometric_quantities(params['meshfile'])
projector = df.Identity(3) - df.outer(normal, normal)
TFS = df.TensorFunctionSpace(mesh, 'P', 1)
p = polarity_functions[-1]
grad_p = df.project(projector * df.grad(p) * projector, TFS)
div_p = df.tr(grad_p)
omega = df.skew(grad_p)
Theta = df.assemble(df.inner(div_p, div_p) * df.dx(mesh))
Phi = df.assemble(df.inner(omega, omega) * df.dx(mesh))
defect_locations, defect_eigvals = find_defects(mesh, p, grad_p, pmag[-1])
assert len(defect_locations)==2
params['Theta_norm'] = Theta
params['Phi_norm'] = Phi
# clean up
for key in ["converged", "dtpnorm_ave"]:
if key in params.keys():
del params[key]
with open(os.path.join(DIR, params['timestamp'] + '_polarity.json'), "w") as fp:
json.dump(params, fp, indent=4)
plotter = vd.plotter.Plotter(axes=4)
poly = vd.utils.buildPolyData(mesh.coordinates(), mesh.cells())
scalar_actor = vd.mesh.Mesh(poly)
scalar_actor.computeNormals(points=True, cells=True)
scalar_actor.cmap(scalar_cmap, pmag[0], vmin=pmin, vmax=pmax, n=n_cmap_vals)
scalar_actor.addScalarBar(
title = r'$\sqrt{\left(\frac{\beta}{\alpha}\right)} \; \left\vert\mathbf{p}\right\vert$',
pos=(0.8, 0.1), nlabels=2,
titleYOffset=15, titleFontSize=28, size=(100, 600))
scalar_actor.scalarbar.SetLabelFormat("%0.1f")
plotter.add(scalar_actor)
defects = vd.shapes.Spheres(defect_locations, c='red', r=0.02)
plotter.add(defects)
vec = defect_locations[1] - defect_locations[0]
uvec = vec / np.linalg.norm(vec)
line = vd.Line(-1.1*uvec, 1.1*uvec, lw=6, c='black')
plotter.add(line)
endPoints = mesh.coordinates() + vector_scale * polarity[0]
vector_actor = vd.shapes.Arrows(mesh.coordinates(), endPoints)
vector_actor.color(vector_color)
plotter.add(vector_actor)
def update(idx):
nonlocal vector_actor
scalar_actor.pointdata['polarity'] = pmag[idx]
scalar_actor.cmap(scalar_cmap, pmag[idx], vmin=pmin, vmax=pmax, n=n_cmap_vals)
plotter.remove(vector_actor)
endPoints = mesh.coordinates() + vector_scale * polarity[idx]
vector_actor = vd.shapes.Arrows(mesh.coordinates(), endPoints)
vector_actor.color(vector_color)
plotter.add(vector_actor)
def slider_update(widget, event):
value = widget.GetRepresentation().GetValue()
idx = (abs(times-value)).argmin()
update(idx)
slider = plotter.addSlider2D(slider_update, pos=[(0.25,0.05), (0.75,0.05)],
xmin=times[0], xmax=times.max(),
value=times[0], title=r"$\alpha \, t$")
slider.GetRepresentation().SetLabelFormat('%0.1f')
# optimized camera position
cam = dict(pos=(4.952, 3.107, 1.804),
focalPoint=(0.07537, 0.1585, -0.07462),
viewup=(-0.2810, -0.1402, 0.9494),
distance=6.000,
clippingRange=(2.836, 10.26))
vd.show(interactive=(not offscreen), camera=cam, zoom=1.4)
# make movie
if offscreen:
print('Saving movie...')
FPS = 10
movFile = os.path.join(DIR, '%s_polarity.mov' % params['timestamp'])
fps = float(FPS)
command = "ffmpeg -y -r"
options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
tmp_dir = TemporaryDirectory()
get_filename = lambda x: os.path.join(tmp_dir.name, x)
for tt in progressbar.progressbar(range(len(times))):
idx = (abs(times-times[tt])).argmin()
update(idx)
slider.GetRepresentation().SetValue(times[tt])
fr = get_filename("%03d.png" % tt)
vd.io.screenshot(fr)
os.system(command + " " + str(fps)
+ " -i " + tmp_dir.name + os.sep
+ "%03d.png " + options + " " + movFile)
tmp_dir.cleanup()
if __name__ == '__main__':
import argparse, json
parser = argparse.ArgumentParser()
parser.add_argument('-j','--jsonfile', help='parameters file', required=True)
parser.add_argument('--onscreen', action=argparse.BooleanOptionalAction)
args = parser.parse_args()
with open(args.jsonfile) as jsonFile:
params = json.load(jsonFile)
(times, geometry, topology, concentration) = get_data(params, DIR='')
visualize(params, DIR=os.path.dirname(args.jsonfile), offscreen=(not args.onscreen))
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment