# PDESolvers in CompuCell3D¶

One of the most important and time consuming parts of the CC3D simulation is to solve all sorts of Partial Differential Equations which describe behavior of certain simulation objects (usually chemical fields). Most of the CC3D PDE solvers solve PDE with diffusive terms. Because we are dealing with moving boundary condition problems it was easiest and probably most practical to use explicit scheme of Finite Difference method. Most of CC3D PDE solvers run on multi core architectures and we also have GPU solvers which run and high performance GPU’s and they also provide biggest speedups in terms of performance. Because CC3D solvers were implemented at different CC3D life cycle and often in response to particular user requests, CC3DML specification may differ from solver to solver. However, the basic structure of CC3DML PDE solver code follows similar pattern .

- FlexibleDiffusionSolver
- Instabilities of the Forward Euler Method
- DiffusionSolverFE
- AdvectionDiffusionSolver.
- FastDiffusionSolver2D
- kernel_diffusion_solver_2D
- ReactionDiffusionSolver
- Steady State diffusion solver

## FlexibleDiffusionSolver¶

This steppable is one of the basic and most important modules in CompuCell3D simulations.

Tip

Starting from version `3.6.2`

we developed `DiffusionSolverFE`

which eliminates several inconveniences of `FlexibleDiffusionSolver`

.

As the name suggests it is responsible for solving diffusion equation but in addition to this it also handles chemical secretion which maybe thought of as being part of general diffusion equation:

where `k`

is a decay constant of concentration `c`

and `D`

is the
diffusion constant. The term called `secretion`

has the meaning as
described below.

Example syntax for FlexibleDiffusionSolverFE looks as follows:

```
<Steppable Type="FlexibleDiffusionSolverFE">
<AutoscaleDiffusion/>
<DiffusionField Name="FGF8">
<DiffusionData>
<FieldName>FGF8</FieldName>
<DiffusionConstant>0.1</DiffusionConstant>
<DecayConstant>0.002</DecayConstant>
<ExtraTimesPerMCS>5</ExtraTimesPerMCS>
<DeltaT>0.1</DeltaT>
<DeltaX>1.0</DeltaX>
<DoNotDiffuseTo>Bacteria</DoNotDiffuseTo>
<InitialConcentrationExpression>x*y
</InitialConcentrationExpression>
</DiffusionData>
<SecretionData>
<Secretion Type="Amoeba">0.1</Secretion>
</SecretionData>
<BoundaryConditions>
<Plane Axis="X">
<ConstantValue PlanePosition="Min" Value="10.0"/>
<ConstantValue PlanePosition="Max" Value="10.0"/>
</Plane>
<Plane Axis="Y">
<ConstantDerivative PlanePosition="Min" Value="10.0"/>
<ConstantDerivative PlanePosition="Max" Value="10.0"/>
</Plane>
</BoundaryConditions>
</DiffusionField>
<DiffusionField Name="FGF">
<DiffusionData>
<FieldName>FGF</FieldName>
<DiffusionConstant>0.02</DiffusionConstant>
<DecayConstant>0.001</DecayConstant>
<DeltaT>0.01</DeltaT>
<DeltaX>0.1</DeltaX>
<DoNotDiffuseTo>Bacteria</DoNotDiffuseTo>
</DiffusionData>
<SecretionData>
<SecretionOnContact Type="Medium"
SecreteOnContactWith="Amoeba">0.1</SecretionOnContact>
<Secretion Type="Amoeba">0.1</Secretion>
</SecretionData>
</DiffusionField>
</Steppable>
```

We define sections that describe a field on which the steppable is to
operate. In our case we declare just two diffusion fields - `FGF8`

and `FGF`

.

Warning

When you want to solve more than one diffusive field using given diffusion
you should declare those multiple fiuelds within a single definition of the diffusion solver. You cannot
include two diffusion solver definitions - one with e.g. `FGF8`

and another with `FGF`

.
In other words you can only define e.g. `FlexibleDiffusionSolverFE`

once per simulation.
You can however use FlexibleDiffusionSolverFE for e.g. `FGF`

and `DiffusionSolverFE`

for e.g. `FGF8`

Inside the diffusion field we specify sections describing diffusion and
secretion. Let’s take a look at `DiffusionData`

section first:

```
<DiffusionField Name="FGF8">
<DiffusionData>
<FieldName>FGF8</FieldName>
<DiffusionConstant>0.1</DiffusionConstant>
<DecayConstant>0.002</DecayConstant>
<ExtraTimesPerMCS>5</ExtraTimesPerMCS>
<DeltaT>0.1</DeltaT>
<DeltaX>1.0</DeltaX>
<DoNotDiffuseTo>Bacteria</DoNotDiffuseTo>
<InitialConcentrationExpression>x*y
</InitialConcentrationExpression>
</DiffusionData>
```

We give a name (`FGF8`

) to the diffusion field – this is required as we
will refer to this field in other modules.

Next we specify diffusion constant and decay constant.

Warning

We use Forward Euler Method to solve these equations. This is not a stable method for solving diffusion equation and we do not perform stability checks. If you enter too high diffusion constant for example you may end up with unstable (wrong) solution. Always test your parameters to make sure you are not in the unstable region.

We may also specify cells that will not participate in the diffusion.
You do it using `<DoNotDiffuseTo>`

tag. In this example we do not let any `FGF`

diffuse
into `Bacteria`

cells. You may of course use as many as necessary
`<DoNotDiffuseTo>`

tags. To prevent decay of a chemical in certain cells
we use syntax:

```
<DoNotDecayIn>Medium</DoNotDecayIn>
```

In addition to diffusion parameters we may specify how secretion should
proceed. `SecretionData`

section contains all the necessary information to
tell CompuCell how to handle secretion. Let’s study the example:

```
<SecretionData>
<SecretionOnContact Type="Medium" SecreteOnContactWith="Amoeba">0.1</SecretionOnContact>
<Secretion Type="Amoeba">0.1</Secretion>
</SecretionData>
```

Here we have a definition two major secretion modes. Line:

```
<Secretion Type="Amoeba">0.1</Secretion>
```

ensures that every cell of type `Amoeba`

will get `0.1`

increase in
concentration every MCS. Line:

```
<SecretionOnContact Type="Medium" SecreteOnContactWith="Amoeba">0.1</SecretionOnContact>
```

means that cells of type `Medium`

will get additional `0.1`

increase in
concentration **but only when** they touch cell of type `Amoeba`

. This mode of
secretion is called `SecretionOnContact`

.

We can also see new CC3DML tags `<DeltaT>`

and `<DeltaX>`

. Their values
determine the correspondence between MCS and actual time and between
lattice spacing and actual spacing size. In this example for the first
diffusion field one MCS corresponds to ```
0.``1 units of actual time and
lattice spacing is equal ``1
```

unit of actual length. What is happening here
is that the diffusion constant gets multiplied by:

```
DeltaT/(DeltaX * DeltaX)
```

provided the decay constant is set to 0. If the decay constant is not zero DeltaT appears additionally in the term (in the explicit numerical approximation of the diffusion equation solution) containing decay constant so in this case it is more than simple diffusion constant rescaling.

`DeltaT`

and `DeltaX`

settings are closely related to `ExtraTimesPerMCS`

setting which allows calling of diffusion (and only diffusion) more than
once per MCS. The number of extra calls per MCS is specified by the user
on a per-field basis using `ExtraTimesPerMCS`

tag.

Warning

When using `ExtraTimesPerMCS`

secretion functions will
called only once per MCS. This is different than using `PDESolverCaller`

where entire module is called multiple times (this include diffusion and
secretion for all fields).

Tip

We recommend that you stay away from redefining `DeltaX`

and
`DeltaT`

and assume that your diffusion/decay coefficients are expressed
in units of pixel (distance) and MCS (time). This way when you assign
physical time and distance units to MCS and pixels you can easily
obtain diffusion and decay constants. `DeltaX`

and `DeltaT`

introduce
unnecessary complications.

The `AutoscaleDiffusion`

tag tells CC3D to automatically rescale diffusion
constant when switching between square and hex lattices. In previous
versions of CC3D such scaling had to be done manually to ensure that
solutions diffusion of equation on different lattices match. Here we
introduced for user convenience a simple tag that does rescaling
automatically. The rescaling factor comes from the fact that the
discretization of the divergence term in the diffusion equation has
factors such as unit lengths, using surface are and pixel/voxel volume
in it. On a square lattice all those values have numerical value of `1.0`

.
On hex lattice, and for that matter of non-square latticeses, only
pixel/voxel volume has numerical value of `1`

. All other quantities have
values different than `1.0`

which causes the necessity to rescale
diffusion constant. The detail of the hex lattice derivation will be
presented in the “Introduction to Hexagonal Lattices in CompuCell3D” -
http://www.compucell3d.org/BinDoc/cc3d_binaries/Manuals/HexagonalLattice.pdf .

The `FlexibleDiffusionSolver`

is also capable of solving simple coupled
diffusion type PDE of the form:

where \(m_c\), \(m_f\), \(n_c\) , \(n_f\), \(p_c\), \(p_d\) are coupling coefficients. To code the above equations in CC3DML syntax you need to use the following syntax:

```
<Steppable Type="FlexibleDiffusionSolverFE">
<DiffusionField>
<DiffusionData>
<FieldName>c</FieldName>
<DiffusionConstant>0.1</DiffusionConstant>
<DecayConstant>0.002</DecayConstant>
<CouplingTerm InteractingFieldName=”d” CouplingCoefficent=”0.1”/>
<CouplingTerm InteractingFieldName=”f” CouplingCoefficent=”0.2”/>
<DeltaT>0.1</DeltaT>
<DeltaX>1.0</DeltaX>
<DoNotDiffuseTo>Bacteria</DoNotDiffuseTo>
</DiffusionData>
<SecretionData>
<Secretion Type="Amoeba">0.1</Secretion>
</SecretionData>
</DiffusionField>
<DiffusionField>
<DiffusionData>
<FieldName>d</FieldName>
<DiffusionConstant>0.02</DiffusionConstant>
<DecayConstant>0.001</DecayConstant>
<CouplingTerm InteractingFieldName=”c” CouplingCoefficent=”-0.1”/>
<CouplingTerm InteractingFieldName=”f” CouplingCoefficent=”-0.2”/>
<DeltaT>0.01</DeltaT>
<DeltaX>0.1</DeltaX>
<DoNotDiffuseTo>Bacteria</DoNotDiffuseTo>
</DiffusionData>
<SecretionData>
<Secretion Type="Amoeba">0.1</Secretion>
</SecretionData>
</DiffusionField>
<DiffusionField>
<DiffusionData>
<FieldName>f</FieldName>
<DiffusionConstant>0.02</DiffusionConstant>
<DecayConstant>0.001</DecayConstant>
<CouplingTerm InteractingFieldName=”c” CouplingCoefficent=”-0.2”/>
<CouplingTerm InteractingFieldName=”d” CouplingCoefficent=”0.2”/>
<DeltaT>0.01</DeltaT>
<DeltaX>0.1</DeltaX>
<DoNotDiffuseTo>Bacteria</DoNotDiffuseTo>
</DiffusionData>
<SecretionData>
<Secretion Type="Amoeba">0.1</Secretion>
</SecretionData>
</DiffusionField>
</Steppable>
```

As one can see the only addition that is required to couple diffusion equations has simple syntax:

```
<CouplingTerm InteractingFieldName=”c” CouplingCoefficent=”-0.1”/>
<CouplingTerm InteractingFieldName=”f” CouplingCoefficent=”-0.2”/>
```

## 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
- If

you will need to to call solver multiple times per MCS.

- Set
`<ExtraTimesPerMCS>`

to N-1 where:

and

## Initial Conditions¶

We can specify initial concentration as a function of `x`

, `y`

, `z`

coordinates using `<InitialConcentrationExpression>`

tag we use `muParser`

syntax to type the expression. Alternatively we may use
`ConcentrationFileName`

tag to specify a text file that contains values of
concentration for every pixel:

```
<ConcentrationFileName>initialConcentration2D.txt</ConcentrationFileName>
```

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
has, say `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*
```

where `x`

, `y`

, `z`

, denote coordinate of the pixel. `c`

is the value of the
concentration.

**Example:**

```
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.

## Boundary Conditions¶

All standard solvers (`Flexible`

, `Fast`

, and `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
inside `<DiffusionField>`

tags.

The `<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
boundary conditions. `PlanePosition='Min"`

denotes plane parallel to `yz`

plane passing through `x=0`

. Similarly `PlanePosition="Min"`

denotes plane
parallel to `yz`

plane passing through `x=fieldDimX-1`

where `fieldDimX`

is `x`

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
at `y=fieldFimY-1`

we have dirichlet boundary condition.

To specify periodic boundary conditions along, say, `x`

axis we use the
following syntax:

```
<Plane Axis="X">
<Periodic/>
</Plane>
```

Notice, that `<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 `x=fieldDimX-1`

.

## DiffusionSolverFE¶

DiffusionSolverFE is new solver in 3.6.2 and is intended to fully
replace `FlexibleDiffusionSolverFE`

. It eliminates several limitations
and inconveniences of `FlexibleDiffusionSolverFE`

and provides new
features such as GPU implementation or cell type dependent
diffusion/decay coefficients. In addition it also eliminates the need
to rescale diffusion/decay/secretion constants. It checks stability
condition of the PDE and then rescales appropriately all coefficients
and computes how many extra times per MCS the solver has to be called. It
makes those extra calls automatically.

Warning

One of the key differences between `FlexibleDiffusionSolverFE`

and
`DiffusionSolverFE`

is the way in which secretion is treated. In
`FlexibleDiffusionSolverFE`

all secretion amount is done once followed by
possibly multiple diffusion calls to diffusion (to avoid numerical
instabilities). In `DiffusionSolverFE`

the default mode of operation is
such that multiple secretion and diffusion calls are interleaved.
This means that instead of secreting full amount for a given MCS and
diffusing it, the `DiffusionSolverFE`

secretes substance gradually so that
there is equal amount of secretion before each call of the diffusion.
One can change this behavior by adding `<DoNotScaleSecretion/>`

to
definition of the diffusion solver e.g.

With such definition the `DiffusionSolverFE`

will behave like
`FlexibleDiffusionSolverFE`

as far as computation.

Note

`DiffusionSolverFE`

autoscales diffusion discretization
depending on the lattice so that `<AutoscaleDiffusion/>`

we used in
`FlexibleDiffusionSolverFE`

is unnecessary.
This may result in slow performance so users have to be aware that those
extra calls to the solver may be the cause.

Typical syntax for the `DiffusionSolverFE`

may look like example below:

```
<Steppable Type="DiffusionSolverFE">
<DiffusionField Name="ATTR">
<DiffusionData>
<FieldName>ATTR</FieldName>
<GlobalDiffusionConstant>0.1</GlobalDiffusionConstant>
<GlobalDecayConstant>5e-05</GlobalDecayConstant>
<DiffusionCoefficient CellType="Red">0.0</DiffusionCoefficient>
</DiffusionData>
<SecretionData>
<Secretion Type="Bacterium">100</Secretion>
</SecretionData>
<BoundaryConditions>
<Plane Axis="X">
<Periodic/>
</Plane>
<Plane Axis="Y">
<Periodic/>
</Plane>
</BoundaryConditions>
</DiffusionField>
</Steppable>
```

The syntax resembles the syntax for `FlexibleDiffusionSolverFE`

. We
specify global diffusion constant by using `<GlobalDiffusionConstant>`

tag. This specifies diffusion coefficient which applies to entire region
of the simulation. We can override this specification for regions
occupied by certain cell types by using the following syntax:

```
<DiffusionCoefficient CellType="Red">0.0</DiffusionCoefficient>
```

Similar principles apply to decay constant and we use
`<GlobalDecayConstant>`

tag to specify global decay coefficient and

to override global definition for regions occupied by Red cells.

We do not support `<DeltaX>`

, `<DeltaT>`

or `<ExtraTimesPerMCS>`

tags.

Note

`DiffusionSolverFE`

autoscales diffusion discretization
depending on the lattice so that `<AutoscaleDiffusion/>`

we used in
`FlexibleDiffusionSolverFE`

is unnecessary.

### Running DiffusionSolver on GPU¶

To run `DiffusionSolverFE`

on GPU all we have to do (besides having OpenCL
compatible GPU and correct drives installed) to replace first line of
solver specification:

```
<Steppable Type="DiffusionSolverFE">
```

with

```
<Steppable Type="DiffusionSolverFE_OpenCL">
```

Note

Depending on your computer hardware you may or may not be able to take advantage of GPU capabilities.

## AdvectionDiffusionSolver.¶

Note

This is an experimental module and was not fully curated.

`AdvectionDiffusionSolver`

solves advection diffusion equation on a cell field as
opposed to grid. Of course, the inaccuracies are bigger than in the case
of PDE being solved on the grid but on the other hand solving the PDE on
a cell field means that we associate concentration with a given cell (not
just with a lattice point). This means that as cells move so does the
concentration. In other words we get advection for free. The
mathematical treatment of this kind of approximation was described in
Phys. Rev. E 72, 041909 (2005) paper by D.Dan et al.

The equation solved by this steppable is of the type:

where \(c\) denotes concentration , \(D\) is diffusion constant, \(k\) decay constant, \(\vec{\nu}\) is velocity field.

In addition to just solving advection-diffusion equation this module
allows users to specify secretion rates of the cells as well as
different secretion modes. The syntax for this module follows same
pattern as `FlexibleDiffusionSolverFE`

.

Example syntax:

```
<Steppable Type="AdvectionDiffusionSolverFE">
<DiffusionField Name="FGF">
<DiffusionData>
<FieldName>FGF</FieldName>
<DiffusionConstant>0.05</DiffusionConstant>
<DecayConstant>0.003</DecayConstant>
<ConcentrationFileName>flowFieldConcentration2D.txt</ConcentrationFileName>
<DoNotDiffuseTo>Wall</DoNotDiffuseTo>
</DiffusionData>
<SecretionData>
<Secretion Type="Fluid">0.5</Secretion>
<SecretionOnContact Type="Fluid"
<SecreteOnContactWith="Wall">0.3</SecretionOnContact>
</SecretionData>
</DiffusionField>
</Steppable>
```

## FastDiffusionSolver2D¶

Note

The current implementation of `DiffusionSolverFE`

may actually be faster than this legacy module. So
before commiting to it please verify the performance of teh `DiffusionSolverFE`

vs `FastDiffusionSolverFE`

`FastDiffusionSolver2DFE`

module is a simplified version of the
`FlexibleDiffusionSolverFE`

steppable. It runs several times faster that
flexible solver but lacks some of its features. Typical syntax is shown
below:

```
<Steppable Type="FastDiffusionSolver2DFE">
<DiffusionField Name="FGF">
<DiffusionData>
<UseBoxWatcher/>
<FieldName>FGF</FieldName>
<DiffusionConstant>0.010</DiffusionConstant>
<DecayConstant>0.003</DecayConstant>
<ExtraTimesPerMCS>2</ExtraTimesPerMCS>
<DoNotDecayIn>Wall</DoNotDecay>
<ConcentrationFileName>Demos/diffusion/diffusion_2D_fast_box.pulse.txt
</ConcentrationFileName>
</DiffusionData>
</DiffusionField>
</Steppable>
```

In particular, for fast solver you cannot specify cells into which diffusion is prohibited. However, you may specify cell types where diffusant decay is prohibited

For explanation on how `<ExtraTimesPerMCS>`

tag works works see section on
`FlexibleDiffusionSolverFE`

.

## KernelDiffusionSolver¶

This diffusion solver has the advantage over previous solvers that it
can handle large diffusion constants. It is also stable. However, it
does not accept options like `<DoNotDiffuseTo>`

or `<DoNotDecayIn>`

. It also
requires periodic boundary conditions.

Simply put KernelDiffusionSolver solves diffusion equation:

with fixed, periodic boundary conditions on the edges of the lattice.
This is different from `FlexibleDiffusionSolverFE`

where the boundary
conditions evolve. You also need to choose a proper `Kernel range (K)`

according to the value of diffusion constant. Usually when \(K^2e^{-K^2/{4D}}\)
is small (this is the main part of the
integrand), the approximation converges to the exact value.

The syntax for this solver is as follows:

```
<Steppable Type="KernelDiffusionSolver">
<DiffusionField Name="FGF">
<Kernel>4</Kernel>
<DiffusionData>
<FieldName>FGF</FieldName>
<DiffusionConstant>1.0</DiffusionConstant>
<DecayConstant>0.000</DecayConstant>
<ConcentrationFileName>
Demos/diffusion/diffusion_2D.pulse.txt
</ConcentrationFileName>
</DiffusionData>
</DiffusionField>
</Steppable>
```

Inside `<DiffusionField>`

tag one may also use option `<CoarseGrainFactor>`

. For example:

```
<Steppable Type="KernelDiffusionSolver">
<DiffusionField Name="FGF">
<Kernel>4</Kernel>
<CoarseGrainFactor>2</CoarseGrainFactor>
<DiffusionData>
<FieldName>FGF</FieldName>
<DiffusionConstant>1.0</DiffusionConstant>
<DecayConstant>0.000</DecayConstant>
<ConcentrationFileName>
Demos/diffusion/diffusion_2D.pulse.txt
</ConcentrationFileName>
</DiffusionData>
</DiffusionField>
</Steppable>
```

## ReactionDiffusionSolver¶

The reaction diffusion solver solves the following system of N reaction diffusion equations:

where `W`

denotes cell type

Let’s consider a simple example of such system:

It can be coded as follows:

```
<Steppable Type="ReactionDiffusionSolverFE">
<AutoscaleDiffusion/>
<DiffusionField Name="F">
<DiffusionData>
<FieldName>F</FieldName>
<DiffusionConstant>0.010</DiffusionConstant>
<ConcentrationFileName>
Demos/diffusion/diffusion_2D.pulse.txt
</ConcentrationFileName>
<AdditionalTerm>-0.01*H</AdditionalTerm>
</DiffusionData>
</DiffusionField>
<DiffusionField Name="H">
<DiffusionData>
<FieldName>H</FieldName>
<DiffusionConstant>0.0</DiffusionConstant>
<AdditionalTerm>0.01*F</AdditionalTerm>
</DiffusionData>
</DiffusionField>
</Steppable>
```

Notice how we implement functions `f`

from the general system of
reaction diffusion equations. We simply use `<AdditionalTerm>`

tag and
there we type arithmetic expression involving field names (tags
`<FieldName>`

). In addition to this we may include in those expression
word `CellType`

. For example:

```
<AdditionalTerm>0.01*F*CellType</AdditionalTerm>
```

This means that function `f`

will depend also on `CellType`

. `CellType`

holds the value of the type of the cell at particular location - `x`

, `y`

, `z`

- of the lattice. The inclusion of the cell type might be useful if you
want to use additional terms which may change depending of the cell
type. Then all you have to do is to either use if statements inside
`<AdditionalTerm>`

or form equivalent mathematical expression using
functions allowed by `muParser`

: http://muparser.sourceforge.net/mup_features.html#idDef2

For example, let’s assume that additional term for second equation is the following:

In such a case additional term would be coded as follows:

```
<AdditionalTerm>CellType==1 ? 0.01*F : 0.15*F</AdditionalTerm>
```

Notice that we have used here, so called ternary operator which might be familiar to you from other programing languages such as C or C++ and is equivalent to`` if-then-els``e statement

The syntax of the ternary (aka `if-then-else`

statement) is as follows:

```
condition ? expression if condition is true : expression if condition false
```

Warning

**Important:** If change the above expression to

we will get an XML parsing error. Why? This i because XML parser will think
that `<1`

is the beginning of the new XML element. To fix this you could
use two approaches:

1.Present your expression as `CDATA`

```
<AdditionalTerm>
<![CDATA[
CellType<1 ? 0.01*F : 0.15*F
]]>
</AdditionalTerm>
```

In this case XML parser will correctly interpret the expression enclosed
between `<![CDATA[`

and `]]>`

.

2. Replace XML using equivalent Python syntax - see (http://pythonscriptingmanual.readthedocs.io/en/latest/replacing_cc3dml_with_equivalent_python_syntax.html) in which case you would code the above XML element as the following Python statement:

```
DiffusionDataElmnt\_2.ElementCC3D('AdditionalTerm', {}, 'CellType<1 ? 0.01*F : 0.15*F')
```

The moral from this story is that if like to use muParser in the XML file make sure to use this general syntax:

```
<AdditionalTerm>
<![CDATA[
YOUR EXPRESSION
]]>
</AdditionalTerm>
```

One thing to remember is that computing time of the additional term depends on the level of complexity of this term. Thus, you might get some performance degradation for very complex expressions coded in muParser

Similarly as in the case of `FlexibleDiffusionSolverFE`

we may use
`<AutoscaleDiffusion>`

tag tells CC3D to automatically rescale diffusion
constant. See section `FlexibleDiffusionSolver`

or the `Appendix`

for more
information.

## Steady State diffusion solver¶

Often in the multi-scale simulations we have to deal with chemicals
which have drastically different diffusion constants. For slow diffusion
fields we can use standard explicit solvers (*e.g.* `FlexibleDiffusionSolverFE`

)
but once the diffusion constant becomes large
the number of extra calls to explicit solvers becomes so large that
solving diffusion equation using Forward-Euler based solvers is simply
impractical. In situations where the diffusion constant is so large that
the solution of the diffusion equation is not that much different from
the asymptotic solution (i.e. at \(t=\infty\)) it is often more convenient to use
`SteadyStateDiffusion`

solver which solves Helmholtz equation:

where \(F\) is a source function of the coordinates - it is an input to the equation, \(k\) is decay constant and \(c\) is the concentration. The \(F\) function in CC3D is either given implicitly by specifying cellular secretion or explicitly by specifying concentration \(c\) before solving Helmholtz equation.

The CC3D stead state diffusion solvers are stable and allow solutions for large values of diffusion constants. The example syntax for the steady-state solver is shown below:

```
<Steppable Type="SteadyStateDiffusionSolver2D">
<DiffusionField Name="INIT">
<DiffusionData>
<FieldName>INIT</FieldName>
<DiffusionConstant>1.0</DiffusionConstant>
<DecayConstant>0.01</DecayConstant>
</DiffusionData>
<SecretionData>
<Secretion Type="Body1">1.0</Secretion>
</SecretionData>
<BoundaryConditions>
<Plane Axis="X">
<ConstantValue PlanePosition="Min" Value="10.0"/>
<ConstantValue PlanePosition="Max" Value="5.0"/>
</Plane>
<Plane Axis="Y">
<ConstantDerivaive PlanePosition="Min" Value="0.0"/>
<ConstantDerivaive PlanePosition="Max" Value="0.0"/>
</Plane>
</BoundaryConditions>
</DiffusionField>
</Steppable>
```

The syntax is is similar (actually, almost identical) to the syntax of
the `FlexibleDiffusionSolverFE`

. The only difference is that while
`FlexibleDiffusionSolverFE`

works in in both 2D and 3D users need to
specify the dimensionality of the steady state solver. We use

```
<Steppable Type="SteadyStateDiffusionSolver2D">
```

for 2D simulations when all the cells lie in the `xy`

plane and

```
<Steppable Type="SteadyStateDiffusionSolver">
```

for simulations in 3D.

Note

We can use Python to control secretion in the steady state solvers but
it requires a little bit of low level coding. Implementing secretion
in steady state diffusion solver is different from “regular” Forward
Euler solvers. Steady state solver takes secretion rate that is
specified at `t=0`

and returns the solution at `t=∞`

. For a large diffusion
constants we approximate solution to the PDE during one MCS by using
solution at`` t=∞``. However, this means that if at each MCS secretion
changes we have to do three things 1) zero entire field, 2) set
secretion rate 3) solve steady state solver. The reason we need to zero
entire field is because any value left in the field at `mcs=N`

is
interpreted by the solver as a secretion constant at this location at
`mcs=N+1`

. **Moreover, the the secretion constant needs to have negative
value if we want to secrete positive amount of substance - this weird
requirements comes from the fact that we re using 3:sup:`rd` party
solver which inverts signs of the secretion constants.**

An example below demonstrates how we control secretion of the steady
state in Python. First we need to include tag `<ManageSecretionInPython/>`

in the XML definition of the solver:

```
<Steppable Type="SteadyStateDiffusionSolver2D">
<DiffusionField>
<ManageSecretionInPython/>
<DiffusionData>
<FieldName>FGF</FieldName>
<DiffusionConstant>1.00</DiffusionConstant>
<DecayConstant>0.00001</DecayConstant>
</DiffusionData>
</DiffusionField>
</Steppable>
```

In Python the code to control the secretion involves iteration over
every pixel and adjusting concentration (which as we mentioned will be
interpreted by the solver as a secretion constant) and we have to make
sure that we inherit from `SecretionBasePy`

not `SteppableBasePy`

to ensure
proper ordering of calls to Python module and the C++ diffusion solver.

Note

Make sure you inherit from `SecretionBasePy`

when you try
to manage secretion in the steady state solver using Python. This will
ensure proper ordering of calls to steppable and to C++ solver code.

Note

Once you use `<ManageSecretionInPython/>`

tag in the XML
all secretion tags in the `SecretionData`

will be ignored. In other words,
for this solver you cannot mix secretion specification in Python and
secretion specification in the XML.

```
def __init__(self, _simulator, _frequency=1):
SecretionBasePy.__init__(self, _simulator, _frequency)
def start(self):
self.field = CompuCell.getConcentrationField \
(self.simulator, "FGF")
secrConst = 10
for x, y, z in self.everyPixel(1, 1, 1):
cell = self.cellField[x, y, z]
if cell and cell.type == 1:
self.field[x, y, z] = -secrConst
else:
self.field[x, y, z] = 0.0
def step(self, mcs):
secrConst = mcs
for x, y, z in self.everyPixel(1, 1, 1):
cell = self.cellField[x, y, z]
if cell and cell.type == 1:
self.field[x, y, z] = -secrConst
else:
self.field[x, y, z] = 0.0
```

Warning

Notice that all the pixels that do not secrete **have to be 0.0** as
mentioned above. **If you don’t initialize field values in the
non-secreting pixels to 0.0 you will get wrong results**. The above
code, with comments, is available in our Demo suite (`Demos/SteppableDemos/SecretionSteadyState`

or
`Demos/SteppableDemos/SteadyStateDiffusionSolver`

).