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 On Nov 10, 2020, at 11:33 AM, Henri Colaux <[email protected]<mailto:[email protected]>> wrote: 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). Le mardi 10 novembre 2020 à 16:24:40 UTC+1, [email protected]<http://ens-rennes.fr> a écrit : 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... On 10/11/2020 15:45, 'Guyer, Jonathan E. Dr. (Fed)' via fipy wrote: 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 On Nov 3, 2020, at 12:49 PM, Henri Colaux <[email protected]> wrote: 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 : [X] 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 -- To unsubscribe from this group, send email to [email protected] View this message at https://list.nist.gov/fipy --- To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. <spindiff_v3_6_1_mwe.py> -- To unsubscribe from this group, send email to [email protected] View this message at https://list.nist.gov/fipy --- To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. -- To unsubscribe from this group, send email to [email protected]<mailto:[email protected]> View this message at https://list.nist.gov/fipy <spindiff_v3_6_1_mwe_alt.py> [cid:[email protected]] -- To unsubscribe from this group, send email to [email protected] View this message at https://list.nist.gov/fipy --- To unsubscribe from this group and stop receiving emails from it, send an email to [email protected].
