SteadyStateDiffusionSolver2D Equation

8 months ago by
Hello everyone!

I am working on a project using CompuCell3D, and I am using the PDE Solver "SteadyStateDiffusionSolver2D" to manage the diffusion processes of Oxygen.

However, I can't make sense of the results in the simulation. I have tried to solve the equation given in the documentation (Helmholtz equation) numerically, using Matlab, and the solution is very different.
For instance, in CompuCell3D, there are negative values, which I find unreasonable since I have no cells consuming the substance (there is a single cell diffusing oxygen in the centre of the simulation lattice).
Moreover, I couldn't find any information online regarding the use of the Helmholtz equation (in the form provided in the documentation) for diffusion processes.

Could you please clarify me about this issue or provide me with any documentation that details the methods used in this solver?

Thank you for your time!
This from the manual: Moreover 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 rd party solver which inverts signs of the secretion constants.

Does inverting the sign for the secretion solve the problem for negative values?
written 8 months ago by priyomadhyapok 
Thank you for your suggestion! However, I believe that condition does not apply to my case. I am managing secretion in the XML file and, from what I understood from the Reference Manual, that specification refers only to secretion controlled in Python. In fact, the "SteadyStateDiffusionSolver" Demo that manages secretion in the XML file also has positive values for secretory cells.
written 8 months ago by Francisco Marques 

3 Answers

8 months ago by
Could you post a  simplified version of your simulation?
8 months ago by
This is the code I am using in my simplified version of the simulation (with only one cell in the centre of the lattice):
  <Dimensions x="500" y="300" z="1"/>

<Plugin Name="CellType">
  <CellType Freeze="" TypeId="0" TypeName="Medium"/>
  <CellType Freeze="" TypeId="1" TypeName="Produtor"/>

<Steppable Type="SteadyStateDiffusionSolver2D">
  <!-- Specification of PDE solvers -->
  <DiffusionField Name="O2">
        <Secretion Type="Produtor">1</Secretion>
        <Plane Axis="Y">
           <ConstantDerivative PlanePosition="Max" Value="0.0"/>
           <ConstantDerivative PlanePosition="Min" Value="0.0"/>
        <Plane Axis="X">

(Besides this portion, there's the PIF Initializer)

I am running tests with the Decay Constant as 0 or as 0.39.
When it is is 0, this is the result: Max = 0.0260 and Min = -0.0038

On the other hand, in the Matlab simulation (solving the equation with finite differences and applying Jacobi Method), this is the result: Max = 0.0258 and Min = 0

There are no negative values, although there is a small hole in the centre.

In the simulation with the Decay Constant as 0.39, both graphs have only positive values, but the shapes follow the same pattern as above (peak in CC3D vs cone in Matlab simulation) and the maximum values still differ.

My goal is to understand how the values entered in the XML file evolve during the simulation so that I can determine a realistic set of constants for my actual project.

Thank you!

8 months ago by
When you run Matlab code are you solving Helmholtz equation or are you solving the diffusion equation and run it for a long time to determine the steady state of the solution?

File attached: (7.37 KB)

I am attaching mine example that is very similar to yours. I noticed that that when decay constant is 0 the solver does get unstable. From what I remember with steady state solvers the decay constant had to be different than 0.0. So when I used decay constant of 0.4 the solution looks reasonable.

As a matter of fact the comments in the code for steady state solver state that when decay constant is 0.0 the solution may not exist:

/* elmbda */
/* the constant lambda in the helmholtz equation. if */
/* lambda .gt. 0, a solution may not exist. however, hwscrt will */
/* attempt to find a solution. */

I am also attaching original code of th steady state solver for 2D case. Take a look at the comment section in the biginning portion of the code (line 86 and on)

File attached: HWSCRT.c (81.7 KB)

So in summary , SteadyState solver is an approximation of the diffusion problem for large diffusion constants. There might be some numerical issues when lambda is very small compared to diffusion constants. As far as comparison with Matlab we have not done this. We benchmarked the code agains FEMLAB.

We would love to add more sophisticated PDE solver quite to CC3D at some point because Steady-State solvers are just an approximations and I wonder if the discrepancy you see between CC#D code and Matlab is because of the slightly different formulation of the problem (i.e. diffusion equation vs Helmholtz equation)

Let me know if this is helpful

Thank you for your reply! It helped me to understand the problem better!

I see the difference between those two possibilities. I believe my tests rely on the second option: I am solving the equation D(d2u/dx2 + d2u/dy2) - k * u = F(x,y) for a long time to determine the steady state of the solution.
(it is the equation in the Reference Manual but with the Diffusion Constant in the Laplacian, as it did not converge otherwise).
I understand that this may be a cause for the inconsistent result.

Still, I have two questions:
- In the code of the original solver, it states that "if lambda greater than 0, a solution may not exist", right? (from what I found, ".gt." stands for "greater than") So, shouldn't it still work with lambda = 0? (isn't it Poisson's Equation?)
If lambda <= 0, I understand the equation in the CC3D Reference Manual:(d2u/dx2 + d2u/dy2) - k * u = F(x,y). The minus sign confused me because in every other source there was a plus sign. Now I realise that diffusion processes require a negative sign in the equation.

- Also, there's a difference in the lambda parameter. In many sources, it is squared (for instance 1, 2). However, as seen above, the k parameter in the equation of the Reference Manual (and in the solver comments) is not squared.
That is to say, regarding the meaning of this constant, that when I use 0.39 in my simulation, is it the square of the decay constant that I intend to use? And, therefore, I should input 0.39^2 instead, as it expects the squared value?

Thank you!
written 8 months ago by Francisco Marques 
Please login to add an answer/comment or follow this question.

Similar posts:
Search »