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 Le jeudi 12 novembre 2020 à 19:28:08 UTC+1, jonatha...@nist.gov a écrit : > 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 <colau...@gmail.com> 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, martinu...@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 <colau...@gmail.com> 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 : >> >> >> 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 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> >> >> >> -- >> 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. >> >> >> >> > -- > To unsubscribe from this group, send email to fipy+uns...@list.nist.gov > > View this message at https://list.nist.gov/fipy > > <spindiff_v3_6_1_mwe_alt.py> > > > > -- To unsubscribe from this group, send email to fipy+unsubscr...@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+unsubscr...@list.nist.gov.
