added b field simulation to probe magnet configurations

as5013-test/sim/CMD

+python bfield_fd.py --mode single_magnet -nt 40000 -mu 100 -nx 200 -ny 200
+python bfield_fd.py --mode double_magnet -nt 40000 -mu 100 -nx 200 -ny 200

as5013-test/sim/bfield_fd.py

+#! /usr/bin/env python
+from __future__ import division
+from numpy import *
+import argparse
+from matplotlib import pyplot as plt
+from matplotlib import animation
+plt.style.use('bmh')
+from matplotlib.patches import Rectangle
+import matplotlib.colors as colors
+
+#using sor to solve for electric fields
+eps_0 = 8.85e-12 #permittivity of free space
+mu_0 = 1.25663706e-6 #permittivity of free space (m*kg/s^2/A^2)
+
+fn_ind = 0
+contours = None
+
+def compute_a(mu):
+	#compute square contour integrals for permeability mu
+	a0 = mu[1:,1:] + mu[:-1,1:] + mu[1:,:-1] + mu[:-1,:-1]
+	a1 = .5*(mu[1:,1:] + mu[1:,:-1])
+	a2 = .5*(mu[1:,1:] + mu[:-1,1:])
+	a3 = .5*(mu[:-1,:-1] + mu[:-1,1:])
+	a4 = .5*(mu[:-1,:-1] + mu[1:,:-1])
+	return a0,a1,a2,a3,a4
+def solve_poisson(args,V,rho,a,neuman_bcs=[None,None,None,None]):
+	for ti in range(0,args.nt):
+		V[0] = V[1]
+		if neuman_bcs[0] is not None:
+			V[1,0,:] = V[1,1,:] - neuman_bcs[0]
+		if neuman_bcs[1] is not None:
+			V[1,-1,:] = V[1,-2,:] + neuman_bcs[1]
+		if neuman_bcs[2] is not None:
+			V[1,:,0] = V[1,:,1] - neuman_bcs[2]
+		if neuman_bcs[3] is not None:
+			V[1,:,-1] = V[1,:,-2] + neuman_bcs[3]
+		V[1,1:-1,1:-1] = (1-args.alpha)*V[0,1:-1,1:-1] + \
+				args.alpha*( a[1]*V[0,2:,1:-1] + a[2]*V[0,1:-1,2:] + a[3]*V[0,:-2,1:-1] + a[4]*V[0,1:-1,:-2] )/a[0] + \
+				args.alpha*rho[1:-1,1:-1]/mu_0/a[0]
+
+		R = absolute(V[0,1:-1,1:-1] - V[1,1:-1,1:-1])
+		if amax(R) < 1e-6:
+			print "broke after %d iterations"%ti
+			break
+	return V
+def compute_H(V,dx,dy):
+	#take derivatives of V and shift indices
+	Hx = -.5*(V[-1,1:,1:] - V[-1,1:,:-1])/dx - .5*(V[-1,:-1,1:] - V[-1,:-1,:-1])/dx
+	Hy = -.5*(V[-1,1:,1:] - V[-1,:-1,1:])/dy - .5*(V[-1,1:,:-1] - V[-1,:-1,:-1])/dy
+	Hmag = sqrt(Hx*Hx + Hy*Hy)
+	return Hx,Hy,Hmag
+
+
+def run_single_magnet(args):
+	nx = args.nx; dx = args.xs/nx
+	ny = args.ny; dy = args.ys/ny	
+	py1 = int(args.ny/2 - int(args.l/2/dy)) #index for plate location
+	py2 = int(args.ny/2 + int(args.l/2/dy)) #index for plate location
+	py_sep1 = int(args.ny/2 - int(args.l/2/dy) - int(args.sep/dy))
+	px1 = int(args.nx/2 - int(args.w/2/dx)) #index for plate edge
+	px2 = int(args.nx/2 + int(args.w/2/dx))
+	V = zeros((2,ny,nx)) #use twice time steps for half-stepping
+	V[:,:,0] = zeros(ny)
+	V[:,0,:] = zeros(nx)
+	V[:,:,-1] = zeros(ny)
+	V[:,-1,:] = zeros(nx)
+	
+	mu = ones((ny-1,nx-1))
+	mu[ py1:py2, px1:px2+1 ] = args.mu
+	#mu[ 0:py1, px1:px2+1 ] = args.mu
+
+	divM = zeros_like(V[0]) #divergence of magnetic dipole moment density M = B_r / mu_0
+	vol = args.w*args.w*args.l
+	divM[py1,px1:px2+1] = -args.Br/mu_0/(px2-px1)*vol
+	divM[py2,px1:px2+1] =  args.Br/mu_0/(px2-px1)*vol
+
+	V = solve_poisson(args,V,divM,compute_a(mu),neuman_bcs=[0,0,0,0])
+	Hx,Hy,Hmag = compute_H(V,dx,dy)
+
+
+	fig, axs = plt.subplots(2, 1, gridspec_kw = {'height_ratios':[2.5, 1]})
+	axs[0].set_aspect(1.)
+	x = linspace(-.5*args.xs, .5*args.xs, nx); y = linspace(-.5*args.ys, .5*args.ys, ny);
+	CS = axs[0].contourf(x[1:], y[1:], Hmag, cmag = plt.cm.autumn)
+	#cbar = fig.colorbar(CS, ax=axs[0])
+	#cbar.ax.set_ylabel('field strength A/m')	
+	#plt.contourf(x, y, V[-1])
+	X,Y = meshgrid(x[1:],y[1:])
+	axs[0].set_xlim([x[0],x[-1]])
+	axs[0].set_ylim([args.pf*y[0],args.pf*y[-1]])
+	strm = axs[0].streamplot(X,Y, Hx,Hy,color=Hmag, minlength=args.w/2,linewidth=1, cmap=plt.cm.summer)
+	axs[0].add_patch(Rectangle( (-args.w/2, -args.l/2), args.w, args.l, facecolor="grey", alpha=.4 ))
+	axs[0].plot([-args.xs/2,args.xs/2],[-args.l/2-args.sep,-args.l/2-args.sep],c='k',lw=2,ls='--')
+
+	axs[1].plot(x[:-1],Hy[py_sep1])
+	axs[1].set_xlabel('x (m)')
+	axs[1].set_ylabel('B_z')
+	axs[1].set_yticks([0.])
+	axs[1].set_xlim([-2*args.w,2*args.w]) #set x range based on magnet width
+
+	plt.savefig('run_single_magnet.png', bbox_inches='tight', pad_inches=0)
+	plt.show()
+
+
+def run_double_magnet(args):
+	nx = args.nx; dx = args.xs/nx
+	ny = args.ny; dy = args.ys/ny	
+	py1 = int(args.ny/2 - int(args.l/2/dy)) #index for plate location
+	py2 = int(args.ny/2 + int(args.l/2/dy)) #index for plate location
+	py_sep1 = int(args.ny/2 - int(args.l/2/dy) - int(args.sep/dy))
+	px1 = int(args.nx/2 - int(args.w/dx)) #index for plate edge
+	pxm = int(args.nx/2)
+	px2 = int(args.nx/2 + int(args.w/dx))
+	V = zeros((2,ny,nx)) #use twice time steps for half-stepping
+	V[:,:,0] = zeros(ny)
+	V[:,0,:] = zeros(nx)
+	V[:,:,-1] = zeros(ny)
+	V[:,-1,:] = zeros(nx)
+	
+	mu = ones((ny-1,nx-1))
+	mu[ py1:py2, px1:px2+1 ] = args.mu
+	#mu[ 0:py1, px1:px2+1 ] = args.mu
+
+	divM = zeros_like(V[0]) #divergence of magnetic dipole moment density M = B_r / mu_0
+	vol = args.w*args.w*args.l
+	divM[py1,px1:pxm] = -args.Br/mu_0/(px2-px1)*vol
+	divM[py1,pxm:px2] =  args.Br/mu_0/(px2-px1)*vol
+	divM[py2,px1:pxm] =  args.Br/mu_0/(px2-px1)*vol
+	divM[py2,pxm:px2] = -args.Br/mu_0/(px2-px1)*vol
+
+	V = solve_poisson(args,V,divM,compute_a(mu),neuman_bcs=[0,0,0,0])
+	Hx,Hy,Hmag = compute_H(V,dx,dy)
+
+
+	fig, axs = plt.subplots(2, 1, gridspec_kw = {'height_ratios':[2.5, 1]})
+	axs[0].set_aspect(1.)
+	x = linspace(-.5*args.xs, .5*args.xs, nx); y = linspace(-.5*args.ys, .5*args.ys, ny);
+	CS = axs[0].contourf(x[1:], y[1:], Hmag, cmag = plt.cm.autumn)
+	#cbar = fig.colorbar(CS, ax=axs[0])
+	#cbar.ax.set_ylabel('field strength A/m')	
+	#plt.contourf(x, y, V[-1])
+	X,Y = meshgrid(x[1:],y[1:])
+	axs[0].set_xlim([x[0],x[-1]])
+	axs[0].set_ylim([args.pf*y[0],args.pf*y[-1]])
+	strm = axs[0].streamplot(X,Y, Hx,Hy,color=Hmag, minlength=args.w/2,linewidth=1, cmap=plt.cm.summer)
+	axs[0].add_patch(Rectangle( (-args.w, -args.l/2), .99*args.w, args.l, facecolor="grey", alpha=.4 ))
+	axs[0].add_patch(Rectangle( (0.01*args.w, -args.l/2), args.w, args.l, facecolor="grey", alpha=.4 ))
+	axs[0].plot([-args.xs/2,args.xs/2],[-args.l/2-args.sep,-args.l/2-args.sep],c='k',lw=2,ls='--')
+
+	axs[1].plot(x[:-1],Hy[py_sep1])
+	#axs[1].plot(x[:-1],Hy[py_sep2])
+	#axs[1].plot(x[:-1],Hy[py_sep3])
+	axs[1].set_xlabel('x (m)')
+	axs[1].set_ylabel('B_z')
+	axs[1].set_yticks([0.])
+	axs[1].set_xlim([-2*args.w,2*args.w]) #set x range based on magnet width
+
+	plt.savefig('run_double_magnet.png', bbox_inches='tight', pad_inches=0)
+	plt.show()
+
+
+if __name__ == '__main__':
+	parser = argparse.ArgumentParser()
+	parser.add_argument('-M','--mode',choices=('single_magnet','double_magnet'))
+	parser.add_argument("-nx","--nx",  type=int, default=200, help="grid size")
+	parser.add_argument("-ny","--ny",  type=int, default=200, help="grid size")
+	parser.add_argument("-xs","--xs", type=float, default=.05, help="space extent x")
+	parser.add_argument("-ys","--ys", type=float, default=.04, help="space extent y")
+	parser.add_argument("-pf","--pf", type=float, default=.75, help="plot fraction (how much of simulated space to plot)")
+	parser.add_argument("-l","--l", type=float, default=.01, help="magnet length")
+	parser.add_argument("-w","--w", type=float, default=.002, help="magnet width")
+	parser.add_argument("-sep","--sep", type=float, default=.001, help="separation between magnet bottom and graphing plane")
+	parser.add_argument("-mu","--mu", type=float, default=1., help="relative permeability of magnet")
+	parser.add_argument("-Br","--Br", type=float, default=1., help="residual flux density of magnet (Teslas)")
+	parser.add_argument("-nt","--nt", type=int, default=50, help="num time steps")
+	parser.add_argument("-alpha","--alpha", type=float, default=1., help="alpha, relaxation parameter")
+	args = parser.parse_args()
+	if args.mode == 'single_magnet':
+		run_single_magnet(args)
+	if args.mode == 'double_magnet':
+		run_double_magnet(args)
+	'''elif args.mode == 'moving_dialectric_sweep':
+		Cs = []
+		seps = [(args.ys/args.ny)*i for i in range(3,15)]
+		for sep in seps:
+			args.sep = sep
+			Cs.append(run_moving_dialectric(args))
+			fn_ind += 1
+		print "capacitances: ",list(Cs)
+		print "seps",list(seps)
+		plt.figure()
+		plt.plot(1e3*asarray(seps),1e12*asarray(Cs))
+		plt.title('Adjacent trace Capacitance vs. separation')
+		plt.xlabel('separation distance (mm)')
+		plt.ylabel('capacitance (pF)')
+		plt.ylim([0,1.1*amax(1e12*asarray(Cs))])
+		plt.show()
+	'''
+