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].
# -*- 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
nop = numpy.array([10, 8])
noptt = 50
L = numpy.array([5. , 5., 5.])
dr = .1
dt = .1
echelles_position = dr * numpy.atleast_1d( [ numpy.linspace(0,nop[0],nop[0]) , numpy.linspace(0,nop[1],nop[1]) ])
gr = fipy.Grid2D( nx = nop[0] ,dx = dr, ny = nop[1],dy = dr )
r_simple_fipy = fipy.CellVariable( name="R simple", mesh= gr, value = numpy.tile( echelles_position[1], nop[0] ) )
CLs = []
ordre = {'value' : 0. , 'order' : 1}
CLs.append( fipy.boundaryConditions.NthOrderBoundaryCondition( faces = gr.facesLeft, **ordre) )
CLs.append( fipy.boundaryConditions.NthOrderBoundaryCondition( faces = gr.facesRight, **ordre) )
CLs.append( fipy.boundaryConditions.NthOrderBoundaryCondition( faces = gr.facesTop, **ordre) )
CLs.append( fipy.boundaryConditions.NthOrderBoundaryCondition( faces = gr.facesBottom, **ordre) )
CLs = tuple(CLs)
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 )
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() )
Pol = fipy.CellVariable(name="Aimantation", mesh = gr, hasOld = True, value= 0. )
eqdiff = (fipy.TransientTerm ( r_simple_fipy ) == fipy.DiffusionTerm( r_simple_fipy * diffcoefs ) -\
( Pol - mag_eq ) * (Tau1relaxs * r_simple_fipy ) )
Pol_list = []
for truc in range(noptt):
print(truc)
eqdiff.sweep(var = Pol, boundaryConditions = CLs , dt = dt )
Pol.updateOld()
Pol_list.append( numpy.reshape(Pol, nop ) )
