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:

\begin{eqnarray} \nabla^2c-kc=F \end{eqnarray}

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