Reaction Diffusion System on an Accelerating Domain

In a previous blog post we derived equations
that describe the time behaviour of a reaction-diffusion system on a growing space domain.

In that work we assumed that the velocity of domain growth is an increasing function of
one of the two unknown variables (here, protein concentrations) described by our
reaction diffusion system.

Let us add another layer to this problem and assume that not the growth velocity is
dependend on the unknown protein concentration but that the growth acceleration is
a function of this unknown concentration.

We derived the equations that describe the time dynamics of our previous (velocity) problem
in material coordinates
here
and reproduce them for clarity:

taken with initial and boundary conditions for $u$, $v$, and $g$ as
described earlier.

Let us now introduce a fourth equation to this system to describe the time behaviour of
the growth velocity $a$ as a function of its acceleration:

$$\frac{\partial a}{\partial t} = s(X,t),$$

where $s(X,t)$ is the local growth acceleration.

In this expanded system of partial differential equations the time behaviour of the, now, unknown variable $a$
is described by an uncoupled partial differential equation.
Since there is no spatial coupling in this partial differential equation, we do not need to impose
any boundary conditions.
As initial condition, we will assume that the system is at rest at first:

Notice that we update the ODE stepper of g, a and x_grid every time step with the current state
of their corresponding right-hand-side expressions.
This is done with a call to .set_f_params() on the corresponding objects.

fortiinrange(1,N):sigma_u=sigma_u_func(g)sigma_v=sigma_v_func(g)rho_u=rho_u_func(g,a)rho_v=rho_v_func(g,a)A_u=numpy.diagflat([-sigma_u[j]+rho_u[j]forjinrange(1,J)],-1)+\
numpy.diagflat([1.+sigma_u[0]+rho_u[0]]+[1.+2.*sigma_u[j]forjinrange(1,J-1)]+[1.+sigma_u[J-1]-rho_u[J-1]])+\
numpy.diagflat([-(sigma_u[j]+rho_u[j])forjinrange(0,J-1)],1)B_u=numpy.diagflat([sigma_u[j]-rho_u[j]forjinrange(1,J)],-1)+\
numpy.diagflat([1.-sigma_u[0]-rho_u[0]]+[1.-2.*sigma_u[j]forjinrange(1,J-1)]+[1.-sigma_u[J-1]+rho_u[J-1]])+\
numpy.diagflat([sigma_u[j]+rho_u[j]forjinrange(0,J-1)],1)A_v=numpy.diagflat([-sigma_v[j]+rho_v[j]forjinrange(1,J)],-1)+\
numpy.diagflat([1.+sigma_v[0]+rho_v[0]]+[1.+2.*sigma_v[j]forjinrange(1,J-1)]+[1.+sigma_v[J-1]-rho_v[J-1]])+\
numpy.diagflat([-(sigma_v[j]+rho_v[j])forjinrange(0,J-1)],1)B_v=numpy.diagflat([sigma_v[j]-rho_v[j]forjinrange(1,J)],-1)+\
numpy.diagflat([1.-sigma_v[0]-rho_v[0]]+[1.-2.*sigma_v[j]forjinrange(1,J-1)]+[1.-sigma_v[J-1]+rho_v[J-1]])+\
numpy.diagflat([sigma_v[j]+rho_v[j]forjinrange(0,J-1)],1)U_new=numpy.linalg.solve(A_u,B_u.dot(U)+f_vec_u(U,V,g,a))V_new=numpy.linalg.solve(A_v,B_v.dot(V)+f_vec_v(U,V,g,a))whilea_stepper.successful()anda_stepper.t+dt<ti*dt:a_stepper.integrate(a_stepper.t+dt)whileg_stepper.successful()andg_stepper.t+dt<ti*dt:g_stepper.integrate(g_stepper.t+dt)whilex_stepper.successful()andx_stepper.t+dt<ti*dt:x_stepper.integrate(x_stepper.t+dt)g_stepper=g_stepper.set_f_params(a)x_stepper=x_stepper.set_f_params(a)a_stepper=a_stepper.set_f_params(s(U))U=U_newV=V_new# these are the correct "y" values to save for the current time step since# we integrate only up to t < ti*dtg=g_stepper.ya=a_stepper.yx_grid=x_stepper.yU_record.append(U)V_record.append(V)g_record.append(g)a_record.append(a)x_record.append(x_grid)

As in our
previous code
protein mass is conserved by definition of our boundary conditions for U and V.

To make certain that our numerical integration does conserve protein mass (at least to
some sensible degree) we work out the initial total mass and final total mass respectively: