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!
<Potts> <Dimensions x="500" y="300" z="1"/> <Steps>1000</Steps> <Temperature>10.0</Temperature> <NeighborOrder>2</NeighborOrder> </Potts> <Plugin Name="CellType"> <CellType Freeze="" TypeId="0" TypeName="Medium"/> <CellType Freeze="" TypeId="1" TypeName="Produtor"/> </Plugin> <Steppable Type="SteadyStateDiffusionSolver2D"> <!-- Specification of PDE solvers --> <DiffusionField Name="O2"> <DiffusionData> <FieldName>O2</FieldName> <DiffusionConstant>2500</DiffusionConstant> <DecayConstant>0</DecayConstant> </DiffusionData> <SecretionData> <Secretion Type="Produtor">1</Secretion> </SecretionData> <BoundaryConditions> <Plane Axis="Y"> <ConstantDerivative PlanePosition="Max" Value="0.0"/> <ConstantDerivative PlanePosition="Min" Value="0.0"/> </Plane> <Plane Axis="X"> <Periodic/> </Plane> </BoundaryConditions> </DiffusionField> </Steppable>
(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.
File attached: diffusion_steady_state_ext_potential.zip (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