Instabilities of the Forward Euler Method¶
Most of the PDE solvers in CC3D use Forward Euler explicit numerical scheme. This method is unstable for large diffusion constant. As a matter of fact using D=0.25 with pulse initial condition will lead to instabilities in 2D. To deal with this you would normally use implicit solvers however due to moving boundary conditions that we have to deal with in CC3D simulations, memory requirements, performance and the fact that most diffusion constants encountered in biology are quite low (unfortunately this is not for all chemicals e.g. oxygen ) we decided to use explicit scheme. If you have to use large diffusion constants with explicit solvers you need to do rescaling:
- Set D, Δt, Δx according to your model
you will need to to call solver multiple times per MCS.
<ExtraTimesPerMCS>to N-1 where:
We can specify initial concentration as a function of
<InitialConcentrationExpression> tag we use
syntax to type the expression. Alternatively we may use
ConcentrationFileName tag to specify a text file that contains values of
concentration for every pixel:
The value of concentration of the last pixel read for a given cell
becomes an overall value of concentration for a cell. That is if cell
8 pixels, and you specify different concentration at every
pixel, then cell concentration will be the last one read from the file.
Concentration file format is as follows:
*x y z c*
z , denote coordinate of the pixel.
c is the value of the
0 0 0 1.2 0 0 1 1.4 ...
The initial concentration can also be input from the Python script (typically in the start function of the steppable) but often it is more convenient to type one line of the CC3DML script than few lines in Python.
All standard solvers (
ReactionDiffusion) by default
use the same boundary conditions as the GGH simulation (and those are
specified in the Potts section of the CC3DML script). Users can,
however, override those defaults and use customized boundary conditions
for each field individually. Currently CompuCell3D supports the
following boundary conditions for the diffusing fields: periodic,
constant value (Dirichlet) and constant derivative (von Neumann). To
specify custom boundary condition we include <BoundaryCondition> section
<BoundaryCondition> section describes boundary conditions along
particular axes. For example:
<Plane Axis="X"> <ConstantValue PlanePosition="Min" Value="10.0"/> <ConstantValue PlanePosition="Max" Value="10.0"/> </Plane>
specifies boundary conditions along the
x axis. They are Dirichlet-type
PlanePosition='Min" denotes plane parallel to
plane passing through
PlanePosition="Min" denotes plane
yz plane passing through
dimension of the lattice.
By analogy we specify constant derivative boundary conditions:
<Plane Axis="Y"> <ConstantDerivative PlanePosition="Min" Value="10.0"/> <ConstantDerivative PlanePosition="Max" Value="10.0"/> </Plane>
We can also mix types of boundary conditions along single axis:
<Plane Axis="Y"> <ConstantDerivative PlanePosition="Min" Value="10.0"/> <ConstantValue PlanePosition="Max" Value="0.0"/> </Plane>
Here in the
xz plane at
y=0 we have von Neumann boundary conditions but
y=fieldFimY-1 we have dirichlet boundary condition.
To specify periodic boundary conditions along, say,
x axis we use the
<Plane Axis="X"> <Periodic/> </Plane>
<Periodic> boundary condition specification applies to both
“ends” of the axis i.e. we cannot have periodic boundary conditions at
x=0 and constant derivative at