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] 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 :
>
>
> 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]
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].
# -*- coding: utf-8 -*-
"""
Created on Thu Oct 29 12:44:46 2020
@author: Henri
"""
import numpy
import fipy
def __iter_2Drscale__(rscale):
''' Renvoie un générateur à partir de rscale (en 2D) '''
for idx2, truc1 in enumerate(rscale[1]):
for idx1, truc2 in enumerate(rscale[0]):
yield numpy.array( [ idx1 , idx2 ] ), numpy.array([truc1, truc2])
def __palier_carotte__(rscale,r,V):
tab = numpy.zeros( (len(rscale[0]) , len(rscale[1])), dtype = 'float' )
for idx, truc in __iter_2Drscale__(rscale):
if idx[1] >= r[1]:
tab[idx[0],idx[1]] = V[2]
elif idx[0] < r[0]:
tab[idx[0],idx[1]] = V[0]
else:
tab[idx[0],idx[1]] = V[1]
return tab
# Parameters
nop = numpy.array([10, 8])
noptt = 50
L = numpy.array([5. , 5., 5.])
dr = .1
dt = .01
echelles_position = dr * numpy.atleast_1d( [ numpy.linspace(0,nop[0],nop[0]) , numpy.linspace(0,nop[1],nop[1]) ])
gr = fipy.CylindricalGrid2D( nx = nop[0] ,dx = dr, ny = nop[1],dy = dr )
# Cell variables
diffcoefs = fipy.CellVariable(name = "Diff coef", mesh = gr )
Tau1relaxs = fipy.CellVariable(name = "Relaxation rate", mesh = gr )
mag_eq = fipy.CellVariable(name = "Equilibrium mag",mesh = gr )
# Systeme parameters
diffcoefs.setValue( __palier_carotte__(echelles_position, L, numpy.array([ 1.e-3 , 1.e-3, 1.e-4 ]) ) .ravel() )
Tau1relaxs.setValue( __palier_carotte__(echelles_position, L, numpy.array([ 2. , 2. , 10. ]) ).ravel() )
mag_eq.setValue( __palier_carotte__(echelles_position, L, numpy.array([ 1. , 1. , 10. ]) ).ravel() )
# Équation differentielle
Pol = fipy.CellVariable(name="Aimantation", mesh = gr, hasOld = True, value= 0. )
eqdiff = (fipy.TransientTerm ( ) == fipy.DiffusionTerm( diffcoefs ) -\
( Pol - mag_eq ) * (Tau1relaxs ) )
Pol_list = []
# Sweep
for truc in range(noptt):
print(truc)
eqdiff.sweep(var = Pol, dt = dt )
Pol.updateOld()
Pol_list.append( numpy.reshape(Pol, nop ) )
print(Pol_list[-1])
