<jupyter_start><jupyter_text>Comparing the stray field of a sphere with a dipolar field **Author:** <NAME>**Date:** Jun 2014**Purpose:** We compute the stray field of a uniformly magnetised sphere (using the airbox method) and compare it with the field created by a point dipole with the same magnetic moment. This allows to check how accurate the demag calculation is, and how quickly it converges for finer mesh discretizations.<jupyter_code>from __future__ import division
from finmag.example import sphere_inside_airbox
from finmag.util import helpers
from finmag.util.meshes import mesh_volume, mesh_info, plot_mesh_with_paraview
from finmag.energies.zeeman import DipolarField
import numpy as np
import sys
import finmag
finmag.set_logging_level('INFO')
#np.set_printoptions(precision=1e-8)<jupyter_output>[2014-06-09 18:36:56] INFO: Finmag logging output will be appended to file: '/home/albert/.finmag/global.log'
[2014-06-09 18:36:56] DEBUG: Building modules in 'native'...
[2014-06-09 18:36:57] DEBUG: FinMag          5047:e8225c1b7a79ea431efa470d26532258c63bb6ef
[2014-06-09 18:36:57] DEBUG: Dolfin          1.4.0                Matplotlib      1.3.1               
[2014-06-09 18:36:57] DEBUG: Numpy           1.8.1                Scipy           0.12.0              
[2014-06-09 18:36:57] DEBUG: IPython         2.1.0                Python          2.7.5+              
[2014-06-09 18:36:57] DEBUG: Paraview        4.0.1-1              Sundials        2.5.0               
[2014-06-09 18:36:57] DEBUG: Boost-Python    <unknown>            Linux           Linux Mint 16 Petra 
[2014-06-09 18:36:57] DEBUG: Registering debug signal handler. Press Ctrl-Z at any time to stop execution and jump into the debugger.<jupyter_text>Creating the simulation We begin by defining the relevant parameters. The mesh is a sphere contained inside a box (representing the 'air'), and we define the relevant geometry parameters and magnetisation of the sphere.<jupyter_code>r_sphere = 3.0
r_shell = 30.0
l_box = 100.0
maxh_sphere = 4.0
Ms = 8.6e5
demag_atol = 1e-8
demag_rtol = 1e-8
m_init = [1, 0, 0]
#m_init = [7, -4, 3]  # some random magnetisation orientation

demag_solver_params = {'phi_1': {'absolute_tolerance': demag_atol, 'relative_tolerance': demag_rtol},
                       'phi_2': {'absolute_tolerance': demag_atol, 'relative_tolerance': demag_rtol}}<jupyter_output><empty_output><jupyter_text>Next we create a simulation object with this mesh and the parameters defined above.<jupyter_code>sim = sphere_inside_airbox(r_sphere=r_sphere, r_shell=r_shell, l_box=l_box,
                           maxh_sphere=maxh_sphere, maxh_box=10.0, maxh_shell=10.0,
                           Ms=Ms, m_init=m_init, demag_solver_params=demag_solver_params)<jupyter_output>[2014-06-09 18:37:03] INFO: Finmag logging output will be written to file: '/home/albert/work/code/finmag/doc/ipython_notebooks_src/sphere_inside_airbox.log' (any old content will be overwritten).
[2014-06-09 18:37:03] INFO: Creating Sim object 'sphere_inside_airbox' (rank=0/1).
[2014-06-09 18:37:03] INFO: <Mesh of topological dimension 3 (tetrahedra) with 2360 vertices and 9483 cells, ordered><jupyter_text>We'd like to compare the stray field of the sphere to the field of a point dipole with an equivalent magnetic moment. Ideally, we could just use the volume of a perfect sphere to compute this magnetic moment. However, the spherical mesh inside the box is only a tetrahedral approximation, so the volumes are different.<jupyter_code>mesh_sphere = sim.get_submesh('sphere')
print "Volume of sphere: {}".format(4./3*pi*r_sphere**3)
print "Volume of spherical mesh: {}".format(mesh_volume(mesh_sphere))<jupyter_output>Volume of sphere: 113.097335529
Volume of spherical mesh: 102.753226719<jupyter_text>This is what the mesh of the sphere looks like.<jupyter_code>plot_mesh_with_paraview(mesh_sphere, rescale=0.5)<jupyter_output><empty_output><jupyter_text>Next we create a field representing a point dipole whose magnetic moment is the same as that of the spherical mesh and add it to the simulation.<jupyter_code>vol = mesh_volume(mesh_sphere)
H_dipole = DipolarField(pos=(0, 0, 0), m=m_init, magnitude=Ms*vol)
sim.add(H_dipole)
#print mesh_info(sim.get_submesh('sphere'))<jupyter_output><empty_output><jupyter_text>Plotting the fields<jupyter_code>def probe_field(sim, field_name, N):
    """
    Probe the field on a regular grid of N x N points in the xz-plane.
    Only the 'air' region of the field is probed, all values outside
    that region are masked out.

    """
    coords = sim.mesh.coordinates()
    xmin = coords[:, 0].min(); zmin = coords[:, 0].min()
    xmax = coords[:, 0].max(); zmax = coords[:, 0].max()
    xvals = np.linspace(xmin, zmax, N)
    zvals = np.linspace(zmin, zmax, N)
    pts = np.array([[[x, 0.0, z] for z in zvals] for x in xvals])
    vals = sim.probe_field(field_name, pts, region='air')
    return pts, vals
def plot_field(sim, field_name, N, type='heatmap', figsize=(4, 4)):
    """
    Plot a cross-section of the field in the xz-plane.
    Supported plot types are: 'heatmap' and 'quiver'.

    """
    pts, vals = probe_field(sim, field_name, N)

    X = pts[:, :, 0]  # x-coordinates of mesh vertices
    Z = pts[:, :, 2]  # z-coordinates of mesh vertices
    U = vals[:, :, 0]  # x-coordinates of demag field vectors
    V = vals[:, :, 2]  # z-coordinates of demag field vectors
    
    fig = plt.figure(figsize=figsize)
    ax = fig.gca()
    if type == 'heatmap':
        im = ax.imshow(vals[:, :, 0], cmap=cm.hot)
        fig.colorbar(im)
    elif type == 'quiver':
        C = vals[:, :, 0]  # colour arrows by x-coordinate of field vectors
        im = quiver(X, Z, U, V, C, pivot='middle', cmap=cm.hot)
        ax.set_aspect('equal')
        fig.colorbar(im)<jupyter_output><empty_output><jupyter_text>Stray field of the sphere First we plot the demag field in the air region (i.e., the stray field of the sphere). We show both a heatmap and a vector field plot.<jupyter_code>plot_field(sim, 'Demag', N=100, type='heatmap')
plot_field(sim, 'Demag', N=20, type='quiver', figsize=(6, 6))<jupyter_output><empty_output><jupyter_text>Dipolar field The field of the point dipole looks almost identical, which is a good sign.<jupyter_code>plot_field(sim, 'DipolarField', N=100, type='heatmap')
plot_field(sim, 'DipolarField', N=20, type='quiver', figsize=(6, 6))<jupyter_output><empty_output><jupyter_text>Visualising the errors Next we compute the errors fields and visualise them by various means.<jupyter_code>pts, vals_demag = probe_field(sim, 'Demag', N=100)
pts, vals_dipole = probe_field(sim, 'DipolarField', N=100)
abserr = np.linalg.norm(vals_demag - vals_dipole, axis=2)
relerr = np.ma.masked_invalid(abserr / np.linalg.norm(vals_dipole, axis=2))<jupyter_output><empty_output><jupyter_text>First we show a histogram of the absolute errors. Most of them are very small, which represents the fact that the dipolar field rapidly decays with increasing distance from the sphere.<jupyter_code># IPYTHON_TEST_IGNORE_OUTPUT
hist(abserr.ravel(), bins=30)<jupyter_output><empty_output><jupyter_text>The relative errors are more evenly distributed around a non-zero value.<jupyter_code># IPYTHON_TEST_IGNORE_OUTPUT
hist(relerr[numpy.isfinite(relerr)].ravel(), bins=30)
# IPYTHON_TEST_IGNORE_OUTPUT
print "Maximum absolute difference: {}".format(abserr.max())
print "Mean absolute difference: {}".format(abserr.mean())
print "Maximum relative difference: {}".format(relerr.max())
print "Mean relative difference: {}".format(relerr.mean())<jupyter_output>Maximum absolute difference: 111.113820922
Mean absolute difference: 10.7848269452
Maximum relative difference: 0.494044427848
Mean relative difference: 0.112298876826<jupyter_text>Here is a heatmap of the absolute differences. We can see that the largest differences occur near the sphere itself, as is expected because these are also the locations with the larges field values.<jupyter_code>imshow(abserr, cmap=cm.hot)
colorbar()<jupyter_output><empty_output><jupyter_text>The relative differences show a more interesting distribution. Oddly, there seems to be a small region inside the airbox where the relative differences are virtually zero.<jupyter_code>imshow(relerr, cmap=cm.hot)
colorbar()<jupyter_output><empty_output><jupyter_text>Plotting the errors as a function of mesh discretization Here we compute the absolute and relative errors on a range of increasingly fine meshes. This allows us to see how the errors decay as `maxh` gets smaller.<jupyter_code>def compute_field_err(maxh, atol=1e-6, rtol=1e-6):
    demag_solver_params = {'phi_1': {'absolute_tolerance': atol, 'relative_tolerance': rtol},
                           'phi_2': {'absolute_tolerance': atol, 'relative_tolerance': rtol}}
    sim = sphere_inside_airbox(r_sphere=r_sphere, r_shell=r_shell, l_box=l_box,
                               maxh_sphere=maxh, maxh_box=10.0, maxh_shell=10.0, Ms=Ms,
                               m_init=m_init, demag_solver_params=demag_solver_params)
    print "Sphere volume: {:.2f}".format(mesh_volume(sim.get_submesh('sphere')))
    sim.add(H_dipole)
    vals_demag = sim.get_field_as_dolfin_function('Demag', region='air').vector().array().reshape(3, -1)
    vals_dipole = sim.get_field_as_dolfin_function('DipolarField', region='air').vector().array().reshape(3, -1)
    abserr = np.linalg.norm(vals_demag - vals_dipole, axis=0)
    relerr = abserr / np.linalg.norm(vals_dipole, axis=0)
    return abserr, relerr
abserrs = {}
relerrs = {}
maxhs = [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
for maxh in maxhs:
    if not (maxh in abserrs and maxh in relerrs):
        print("DDD maxh={}".format(maxh))
        sys.stdout.flush()
        abserrs[maxh], relerrs[maxh] = compute_field_err(maxh=maxh, atol=1e-10, rtol=1e-10)
def plot_errs(abs_or_rel, fun):
    errs = {'absolute': abserrs, 'relative': relerrs}[abs_or_rel]
    xvals, yvals = zip(*sorted([(k, fun(v)) for (k, v) in errs.items()]))
    plot(xvals, yvals, 'o-')
    ylim(0, 1.05*max(yvals))
    title('{} of {} error'.format(fun.__name__, abs_or_rel))
    xlabel('maxh')
plot_errs('absolute', max)
plot_errs('absolute', np.mean)
plot_errs('absolute', np.median)
plot_errs('relative', max)
plot_errs('relative', np.mean)
plot_errs('relative', np.median)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11610161305901671, 'fraction_numerical': 0.043465662126919734, 'mean_word_length': 4.350904392764858, 'pattern_counts': {'":': 0, '<': 60, '<?xml version=': 0, '>': 60, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '636956', 'license': 'BSL-1.0', 'n_tokens_mistral': 3908, 'n_tokens_neox': 3529, 'n_words': 997, 'path': 'doc/ipython_notebooks_src/tutorial-comparing-sphere-stray-field-with-dipolar-field.ipynb', 'repo_name': 'davidcortesortuno/finmag', 'size': '10364'}