As part of learning FreeFEN, I tried to model a system of several species ( called f1,f1a,f2,f2a,
fgn, and fb ). The species f1,f2, and fgn are present at some uniform concentration in a liquid
while f1a, f2a, and fb are created at some reaction that occurs in the right i in th e attached
images. Reactions transform f1 → f1a , f2-> f2a, and fgn → fb. The first reaction requires
interaction with a fixed substance on the right, the second requires f1a, and the thrid requires f2a.
The images below show [fgn] being depleted where it can form fb but also with a maxima in a surprising location.
The species diffuse but with a diffusivity that depends on concentration of other species.
It does appear that it is possible for high concntration regions to further concentrate
species due to the diffusivity gradient. Curious if anyone has looked at problems
like this. It does appear that sorting out artifacts from real effects or numerical and math
issues can be quite a challenge 
I sorry, you question a too esoteric, without model, equation, it impossible to say something raisonable.
Thanks. It was just more of an observation than a question. I hesistate to post the
code right now as it is in bad shape and part of the result may be due to bugs.
However, I guess that was my point Seeing an extreme value, in this
case a maximum, in a source free region I though was suggetive or
a math or coding mistake. However, I had a diffusivity dependent on concentration
of other species making some effects like this possible. It is just Fick’s
Law with variable “D” . There exist concentration gradients in two directions
as the reactant fgn is concerted to fb on the right side.
Just curious though if there are similar well studied real systems
that are similar or maybe an example. In the real system the species form
a gel and ultimately a solid that stops diffusion.
My immediate interest is blood coagulation- I have seen some models published
but never looked at the details of hose they did it.
Well, I am not shocked in principle by some density overshoot somewhere.
Without knowing the details, it is hard to say, but naively I would say that it might be a consequence of the boundary conditions: from the pictures I make the following observations
- it looks like outside your droplet there are Dirichlet conditions
- I see some gradients through your interface
Could this overshoot be a consequence of you BC, and some conservation laws playing together?
Thanks. Maybe I’ll set the parameters back to reproduce this and post the code
if there is that much interest. In essence the right side is a wall with a patch
that converts f1 to f1a enzymatically at a non saturating first order rate ( df1a/dt= c* f1 ) .
The other three sides are bulk fluid with constant concentrations of reactants f1,f2, and
fgn with zero f1a,f2a, and fb ( Dirichlet ). The fixed enzyme patch that converts f1 to f1a
I guess has an interface but basically everything is supposed to be continuous
except to the extent the diffusivity approaches zero and creates a solid or gel
I looked at the results in R briefly and the diffusivity gradient was there but total
change was very small only about 2 percent IIRC. I did not numerically check
to see if grad(D) \dot grad(n) compared to D*laplacian(n)
Complete speculation here, but I find your explanation reminiscent of thermal-diffusive instability in flames at non-unity Lewis numbers. Preferential diffusion has some interesting consequences including superadiabatic flame temperature. Such phenomena may be related to what you’re describing.
Thanks, assuming this ref captures the issues you raise,
https://link.springer.com/article/10.1007/s10494-020-00173-7
( and I’m surprised this would ever approach unity lol ),
then there may be similar factors here as I guess diffusivity tends to be theramaly
activated and T gradients to the extent you can define a local temperature,
would produce similar results. I guess doing a microscopic simulation may
be important but right now I’m just looking for math
and coding mistakes
Of course things are more intricate at the microscopic level but you can understand the basic physics using only the kuramoto-sivashinsky equation.
Are these 4th deivatives solvable in FreeFEM?
https://en.wikipedia.org/wiki/Kuramoto%E2%80%93Sivashinsky_equation
and fwiw, I guess a line ratio is a temperature lol.
yes, for bilaplacian use (doc is here)
load “Morley”
for time integration, probably Crank-Nicolson, and forward Richardson extrapolation to treat the nonlinearity…
I guess rather than upload the entire system, I just extracted the part
of interest. This considers the case of diffusion with a non-uniform diffusivity
and diffusivities that produce extreme values in the concentrations in a source
free region. The unknown here is called “f1” and steady state distributions
solved with constant concentration of 100 on each side of a square 2D area.,
This iterates through some diffusivities curves producing varying amounts
of concentration peaks. The first few are kind of “controls” then the
effect is seen more clearly in the latter ones, see the code for details.
fwiw, not sure this is all that profound but I was curious. Maybe someone
can relate it to a well studied issue. I just needed the practice lol.
singd.edp (7.9 KB)
I decided that the easiest way to explore this is just time stepping the equation for
dn/dt but I’m not sure that FreeFem is the fastest way to do that. So, I wrote my own
FD code and got much faster iterations. I changed the geometry slightly and now
have an enzymatic square in the middle of the grid that converts f1 to f1a which
converts f2 to f2a which converts fgn to the gel component fb. I’m not sure how close
this is to steady state but I did get s singificant concentration peak. I want to use the FD results in FreeFEm as initial conidtions but could not find an obvious command to import them. I wrote
a loop ( see the attached edp file ) that sort of reads the txt file and produces the attached
display. However, this is not quite right as the plots from “R” show the square in the center along with the diffusivity plot and concentration of the glue fb. Probably once I get going I will want to have at least two or three methods for sanity check and an FD approach may be one of those
so I expect to do this a bit. What is the best way to load data in this format into a mesh?
fd.edp (368 Bytes)
fddump.zip (545.7 KB)
It looks like triangulate was the way to go here although if you are
looking for features adding column locations for X/y would be nice
The attached medit shot seems to reflect the text file. fwiw, this is a few
thousand time steps of glue or gel deposition around an activating
square. This is a complicated region including a source so its difficult
to say much. This “goo” is used to reduce the diffusivity and it does look
like the FD simulation is a good sanity check on any FEM work. This
wall of “castle” does indeed seem to deplete the precursor inside and prevent
diffusion of the activating enzyme to the outside and just by eye it looks like the
quantities are reasonably conserved- the amount of goo us probably near the amount
of precursor consumed ( transport from the boundaries is probably low ). I initially
forgot the grad(D)grad(N) term and it falls out of the variational form
but seems to be important to getting this result. D only varies may 10 percent or
so but has a steep gradient.
cat fd.edp
load "medit"
// best example I can find I guess,
// https://doc.freefem.org/documentation/mesh-generation.html#the-keyword-%E2%80%9Ctriangulate%E2%80%9D
mesh Th =triangulate("fd2.txt"); // square(sz,sz);
fespace Vh(Th,P1) ;
Vh data;
real[int] mdata(20000);
ifstream ifs("fddump.txt");
//ifs>> data[][0];
data[]=mdata;
real x;
for (int i=0; i<10000; ++i)
{
for (int j=0; j<11; ++j)
{ if (j==8) ifs>> mdata[i];
else ifs>> x;
}
} // i
data[]=mdata;
//cout<<data<<endl;
//plot(data,fill=1,wait=1);
medit("asdf",Th,data);