hi,
I'm a newbie at fipy ; I played with it essentially to solve a poisson equation
of the form "laplacian(a) = j" where grad(a) on the boundary and j are given
by a numpy array (ie essentially non-uniform neumann conditions).
I used to make it work (but maybe not in the proper way) with fipy 2.x
using "FixedFlux" that is now deprecated. I tried to use the faceGrad.constrain
as it seems to be the way to do it now... but after many tries, the poisson
equation is solved in a way that suggest BC are not considered.
The former code is below, as well as the changes made for 3.x version.
If anyone could help, I would really appreciate.
Thank you,
#
#
import numpy as np
import matplotlib.pyplot as plt
from fipy import Grid2D, CellVariable, DiffusionTerm, FixedFlux
#
#
#
def GetBx(B0, L, dl) :
N = [np.int(np.round(i)) for i in np.divide(L, dl)]
Wx = B0*np.linspace(-L[1]/2, L[1]/2, num = N[1]+1)
Bx = np.tile(Wx, (N[0]+1, 1))
return Bx
#
#
def GetBy(B0, L, dl) :
N = [np.int(np.round(i)) for i in np.divide(L, dl)]
Wy = B0*np.linspace(-L[0]/2, L[0]/2, num = N[0]+1)
By = np.repeat(Wy, N[1]+1)
By = np.reshape(By, [N[0]+1, N[1]+1])
return By
#
#
def GetJz(Bx, By, dl) :
dBxdy = np.gradient(Bx, dl[0], dl[1])[1]
dBydx = np.gradient(By, dl[0], dl[1])[0]
Jz = dBydx-dBxdy
return Jz
#
#
def GetAz(Bx, By, Jz, Nc, dl) :
mesh = Grid2D(dx = dl[0], dy = dl[1], nx = Nc[0], ny = Nc[1])
_Az = CellVariable(mesh=mesh, value=0.0)
_Bx = CellVariable(mesh=mesh)
_Bx.value = np.reshape(Bx, Nc[0]*Nc[1], order = 'F')
_By = CellVariable(mesh=mesh)
_By.value = np.reshape(By, Nc[0]*Nc[1], order = 'F')
_Jz = CellVariable(mesh=mesh)
_Jz.value = np.reshape(Jz, Nc[0]*Nc[1], order = 'F')
_Az.equation = (DiffusionTerm(coeff = 1.0)+_Jz == 0)
# beware of the sign of the flux : always consider outward direction
BCs = [FixedFlux(value= _By.getFaceValue(), faces=mesh.getFacesLeft()),
FixedFlux(value=-_By.getFaceValue(), faces=mesh.getFacesRight()),
FixedFlux(value=-_Bx.getFaceValue(), faces=mesh.getFacesBottom()),
FixedFlux(value= _Bx.getFaceValue(), faces=mesh.getFacesTop())]
_Az.equation.solve(var=_Az, boundaryConditions=BCs)
Az = np.reshape(_Az.value, Nc, order = 'F')
return Az
#
#
#
Lb = [20, 10]
dl = [.1, .1]
Nc = [np.int(np.round(i))+1 for i in np.divide(Lb, dl)]
Bx = GetBx(1.0, Lb, dl)
By = GetBy(0.2, Lb, dl)
Jz = GetJz(Bx, By, dl)
Az = GetAz(Bx, By, Jz, Nc, dl)
image = plt.imshow(np.transpose(Az), origin = 'lower')
plt.colorbar(image)
plt.show()
plt.savefig('new.pdf')
To impose the flux on the boundary I tried :
_Az.faceGrad.constrain( _By.arithmeticFaceValue, where = mesh.facesLeft)
_Az.faceGrad.constrain(-_By.arithmeticFaceValue, where = mesh.facesRight)
_Az.faceGrad.constrain(-_Bx.arithmeticFaceValue, where = mesh.facesBottom)
_Az.faceGrad.constrain( _Bx.arithmeticFaceValue, where = mesh.facesTop)
... unsuccessfully.
To understand the physics of the problem, Bx and By are the X & Y components
of the magnetic field, Jz the Z component of the associated current (curl B)
and Az the Z component of the potential vector (B = curl A). With the above
magnetic field topologies, isocontours of Az should be hyperbolas.
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
