Spin diffusion simulation inside a cylinder

[Henri Colaux]

Dear devs,

I am writing a code that simulates the diffusion of magnetisation from proton to proton inside a solid. This process behaves like a standard diffusion program (but with an equilibrium term) so your code have been very helpful and I got some great results so far, but I stambled into some problem whose solution is above my knowledge.

I am trying to simulate a « drill-core » sort of structure with two different components - denoted « domains », one being a small slice inside a much more abundant one, as pictured in the following schematics :

I figured that this system is equivalent to producing a 2D array using Grid2D as picture of the right of this figure, and ensuring that there is no magnetisation transfer (i.e. no flux) taking place at the border of the system. To account for the rotation centre and the cylindrical symmetry, I have taken the effect of the radius directly in the differential equation.Yet, for some reasons, it seems like the boundery conditions is not set correctly in my code, which gave unwanted behaviours.

I have included a minimum working example of the code I am trying to simulate. First and formost, could you let me know if I do everything right for what I am trying to simulate ?

Cheers !

Henri Colaux

[Henri Colaux]

Just in case the figure is not showing :

[Henri Colaux]

One other thing : I tried the « constrain » method with either :

Pol.faceGrad.constrain(((0,),(0,)), where=gr.exteriorFaces)

or

Pol.faceGrad.constrain(((0,),(0,)), where=gr.facesTop)
Pol.faceGrad.constrain(((0,),(0,)), where=gr.facesBottom)
Pol.faceGrad.constrain(((0,),(0,)), where=gr.facesLeft)
Pol.faceGrad.constrain(((0,),(0,)), where=gr.facesRight)

but it does not seem to change anything. Could anyone tell me what I am doing wrong ?
Thank you !

H.

[Guyer, Jonathan E. Dr. (Fed)]

Henri -

My apologies for the delay in answering.

We have an existing mesh generator [`CylindricalGrid2D`](https://www.ctcms.nist.gov/fipy/fipy/generated/fipy.meshes.html#fipy.meshes.factoryMeshes.CylindricalGrid2D) that automatically takes care of calculating cylindrical symmetry. You only need to specify the r and z dimensions of the mesh and FiPy takes care of the wedge-shaped geometry; you do not need to transform your PDE into cylindrical form.

Mathematically, what boundary conditions are you trying to achieve?

I definitely wouldn’t use `NthOrderBoundaryCondition`. Those are only for higher-order PDEs like Cahn-Hilliard/Ginsburg-Landau. Even then, that’s an older way of doing things that we don’t recommend anymore.

`Pol.faceGrad.constrain()` is the syntax we recommend, but you should be aware that as a finite volume code, FiPy has no-flux as the natural boundary condition. While no-flux is not always equivalent to zero normal gradient, for your diffusion equation, they are equivalent. There is no need to apply a constraint if no-flux is what you desire.

- Jon

--
To unsubscribe from this group, send email to fipy+uns...@list.nist.gov
 
View this message at https://list.nist.gov/fipy
---
To unsubscribe from this group and stop receiving emails from it, send an email to fipy+uns...@list.nist.gov.
<spindiff_v3_6_1_mwe.py>

[Martinus WERTS]

Dear Henri,

From your last message I gather that you would like to apply zero-flux conditions on all boundaries, in a cylindrical (r,z) geometry.

In that case, I warmly recommend that you follow Jon's suggestion and use the implicit no-flux boundary condition. Do not explicitly specify a zero-gradient. We had some advection-diffusion test cases here where explicitly specifying zero-flux conditions in Fipy yielded less precise results (loss of mass conservation) than the implicit no-flux.

https://github.com/mhvwerts/MasonWeaver-finite-diff/blob/master/finite-difference-solver-vs-Fipy-finite-volume-solver.pdf


Also, the CylindricalGrid2D is likely exactly what you need for your problem geometry. We used it with success on a system for which an analytic solution exists.

https://github.com/simulkade/FVTool/blob/master/Examples/External/Diffusion1DSpherical_Analytic-vs-FVTool-vs-Fipy/diffusion1Dspherical_analytic_vs_FVTool_vs_Fipy.pdf


Best wishes
Martin


P.S.1. Our example where we used CylindricalGrid2D was created before SphericalGrid1D was implemented in Fipy

P.S.2. For some it is a bit awkward to not specify boundary conditions explicitly. It might be an idea to have a method 'standard_zero_flux_boundaries()' or so, that simply does nothing, but that one could put in one's code so that everything is specified...

[Henri Colaux]

Dear Jon, dear Martinus,

I really appreciate your replies! No worries about the delay. It is a dificult period for all of us :(

Thanks for letting me know that "zero-flux" (the condition I want) is the implicit condition. I have no need to specify it specifically, I just was not sure.

I tried "CylindricalGrid2D" indeed, but I was not sure whether the rotation axes was located at the left of the array, so I decided to use "Grid2D" instead along with a "radius" variable so I am sure of the position of the symmetry axis. Which brings me to my first question: Here is one of the parameter arrays (here, the relaxation time) I use in my MWE, before "ravel()" is apply:

array([[ 2.,  2.,  2.,  2.,  2., 10., 10., 10.],
       [ 2.,  2.,  2.,  2.,  2., 10., 10., 10.],
       [ 2.,  2.,  2.,  2.,  2., 10., 10., 10.],
       [ 2.,  2.,  2.,  2.,  2., 10., 10., 10.],
       [ 2.,  2.,  2.,  2.,  2., 10., 10., 10.],
       [ 2.,  2.,  2.,  2.,  2., 10., 10., 10.],
       [ 2.,  2.,  2.,  2.,  2., 10., 10., 10.],
       [ 2.,  2.,  2.,  2.,  2., 10., 10., 10.],
       [ 2.,  2.,  2.,  2.,  2., 10., 10., 10.],
       [ 2.,  2.,  2.,  2.,  2., 10., 10., 10.]]).ravel()

Could you confirm that the symmetry axis is indeed located on the left of this array ?

In the MWE, I set domains 1 and 2 to have the same caracteristics for simplicity. So, if I was correct of the way my system is programmed, I should have for any time lines with identical values in the polarisation variable "Pol" because there should be a cylindrical symmetry. So it should give identical results to "CylindricalGrid1D". Yet, if I sweep 50 times (for example), here is what I get :

Pol_list[50] =
[[0.79 0.79 0.80 0.99 1.19 9.95 9.95 9.89]
 [0.96 0.80 0.81 0.98 1.19 9.90 9.90 9.89]
 [1.15 0.98 0.81 0.97 1.34 9.91 9.90 9.89]
 [1.14 0.98 0.82 0.98 1.17 9.90 9.90 9.88]
 [1.12 0.98 0.82 0.99 1.17 9.90 9.90 9.95]
 [0.98 0.98 0.82 0.99 1.18 9.90 9.90 9.89]
 [1.15 0.99 0.82 0.99 1.20 9.90 9.90 9.89]
 [1.14 0.98 0.82 0.98 1.34 9.91 9.90 9.89]
 [1.14 0.98 0.81 0.98 1.17 9.90 9.95 9.93]
 [1.13 0.98 0.80 0.80 0.99 9.90 9.95 9.99]]

Here, not only are my lines not identical, but the polarisation goes above the equilibrium value of 1 inside the cylinder - the leftmost values. This does not makes sens unless there is some transfer taking place from the left of the system. Likewise, the rightmost values should by the maximum, yet there are not. I tried with lower time increments and the result stays the same.

I have inclosed the updated MWE with no boundary condition and with CylindricalGrid2D that gave the result above. Could you help my fix it so it gives the same result as would an identical system with CylindricalGrid1D ?

Thanks in advance for your time if you have some :)

PS: I am happy to hear there is now a « SphericalGrid1D » mesh which I missed when starting my program (I have 3.3).

[Martinus WERTS]

Dear Henri,

For setting the values of CellVariables (e.g. initial conditions or variable coefficients), I recommend the following procedure, which avoids thinking about how the finite volume cells are organized in computer memory (and the associated headaches). Try this in a Jupyter Notebook or so.


I assume the follow imports.

import numpy as np
import fipy as fp


Create your cylindrical (r,z) grid. Also, retrieve the 1D arrays which contain the r and z coordinates of each finite-volume cell. You may have noticed that all cells are stored in a 1D array, even for 2D and 3D grids. Note that the coordinates are the coordinates of the cell centers.

gr = fp.CylindricalGrid2D(nr = 10 ,dr = 0.1, nz = 8, dz = 0.1)
rr = gr.cellCenters[0]
zz = gr.cellCenters[1]


Create your relaxation CellVariable. And create a "Viewer" to see what you're doing. (The Viewer was repaired in a recent Fipy release, so please use the latest Fipy)

tau1relaxs = fp.CellVariable(name = "Relaxation rate", mesh = gr )
vw = fp.Viewer(tau1relaxs)


Fill all cells with the value 10.0.

tau1relaxs[:] = 10.


Now the magic begins. There is a syntax for conditionally indexing array elements. The following fills all cells situated at r < 0.5 with the value 2.0. Ask the Viewer object to generate a plot of the CellVariable.

tau1relaxs[rr < 0.5] = 2.
vw.plot()


This can be elaborated further with element-wise 'and' (&).

# make a disc radius 0.7 between z = 0.3 and z = 0.4
tau1relaxs[(rr < 0.7) & (zz > 0.3) & (zz < 0.4)] = 6.
vw.plot()


Hope this helps.

Best regards,
Martin

[Guyer, Jonathan E. Dr. (Fed)]

Henri -

I’m not sure why, as the cells of a CylindricalGrid2D are organized in row-major order, but your coefficients are not being assigned as you intend. See the figure below of Tau1relaxs. My guess is there’s a discrepancy in how you calculate the strides.

I definitely recommend the parametric approach for assigning values advised by Martinus. It’s both less error prone, but will also continue working if you decide to switch to a completely unstructured mesh.

- Jon

<spindiff_v3_6_1_mwe_alt.py>

[Henri Colaux]

Dear Martinus, dear Jon,

Thanks a lot for your replies ! I figured indeed the problem that was highlighted by Jon using the viewer method and was able to fix it by modifing my "__palier_carotte__" function, but it seems a way cleaner way to express my diffusion problem so I am rewritting this part of my code accordingly.

Thank you both for your help ! The only thing I can offer in return though, is only a quotation in our future paper.

I would have one last question if I may: Is there a method that takes the 1D numpy array of CellVariable and returns the 2D or 3D numpy array already in fipy, or must I make one myself ?

Cheers,

Henri Colaux

[Guyer, Jonathan E. Dr. (Fed)]

Henri -

There’s no general facility. We reshape when we need it, e.g., for [`Matplotlib2DGridViewer`](https://github.com/usnistgov/fipy/blob/3cdaccd2b5e963d7d389f5b61b9f5a29ba0762b1/fipy/viewers/matplotlibViewer/matplotlib2DGridViewer.py#L98).

- Jon