Example 1: Demag field in uniformly magnetised sphere

This example computes the demagnetisation field in a uniformly magnetised sphere. We know, of course, that the demag field has to have the opposite direction to the magnetisation, and a magnitude of 1/3 of the magnetisation in this special case.

When using finite element calculations, a crucial (and non-trivial) part of the work is the `finite element mesh generation`_. We provide a very small mesh for this example sphere1.nmesh.h5 which was generated with Netgen_. It describes a sphere of radius 10nm.

Using this, we can write the following nmag script with name sphere1.py:

To execute this script, we have to give its name to the nsim_ executable, for example (on linux):

user@host user> nsim sphere1.py

Let’s discuss the sphere1.py script step by step.

Importing nmag

First we need to import nmag, and any subpackages of nmag that we need (this is just the SI object).

Creating the simulation object

Next, we need to create a simulation object using nmag.Simulation().

Defining (magnetic) materials

After having imported the nmag module into Python’s workspace and after creating the simulation object sim, we need to define a material using nmag.MagMaterial. We are giving it a name (string) which is here "Py" and we are assigning a saturaration magnetisation and an exchange coupling strength. This name of the material chosen here plays an important role as all the `Fields and Subfields`_ will be named taking into account this name. The output files will also use that name to label output data.

All data input and output within nmag happens in SI units by default. To avoid ambiguity here but simultanouesly provide flexibility in the choice of input and output units, an `SI object`_ is provided. This has to be used when defining the material parameters (and in other places). We thus express the saturation magnetisation in Ampere per metre (Ms=SI(1e6,"A/m")) and the exchange coupling constant (often called A in micromagnetism) in Joules per metre (exchange_coupling=SI(13.0e-12, "J/m")).

Loading the mesh

The next step is to load the mesh. The first argument is the file name (sphere1.nmesh.h5). The second argument is a list of tuples. In this example we have a list with one element, and this element is ("sphere", Py). The first part of the tuple ("sphere") is a string (of the user’s choice) and this is the name given to the mesh region 1 (i.e. the space occupied by all simplices that have the region id 1 in the mesh file).

[This information is currently only used for debugging purposes (such as when printing the simulation object).]

The second part of the tuple is the MagMaterial that has been created before and which contains the material properties for this material. In this example, we have assigned the material for PermAlloy to a variable called Py.

[If we had two different materials within the same simulation, then the list would have two entries, i.e. ("sphere", [Dy,Fe2]) in a simulation that takes into account Dy moments and Fe2 moments individually.]

The third argument to load_mesh is an `SI object`_ which defines what real distance should be associated with the length 1.0 as given in the mesh file. In this example, the mesh has been created in nano metres, i.e. the distance 1.0 in the mesh file should correspond to 1 nanometre in the real world. We thus use an SI object of 1nm.

[This defines the internal simulation units.]

Setting the initial magnetisation

To set the initial magnetisation, we use the sim.set_m() method. The field m describes the normalised magnetisation whereas

the field M contains the magnetisation with it’s proper magnitude (i.e. |M| = Ms). For this simulation, we provide a unit vector pointing in x-direction to indicate that we would like the initial magnetisation to point in plus x direction. (We could provide a vector with non-normalised magnitude, which would be normalised automatically. This is convenient for, say, magnetisation pointing 45 degrees between x and y axis: [1,1,0])

Setting the external field

We set the external field using sim.set_H_Ext(). In contrast to sim.set_m(), this methods takes two arguments. The first defines numerical values for the direction and magnitude of the external field. The second determines the meaning of these numerical values using an SI object. Suppose we would like an external field of 1e6A/m acting in the y-direction, then the command would read: sim.set_H_ext([0,1e6,0],SI(1,"A/m")).

Computing the demag field

The sim.advance_time(SI(0,"s")) command is used to compute the demag field, the exchange field, etc and all associated energies, based on the geometry, the magnetisation configuration, and the external field.

If we study dynamical systems, then a call sim.advance_time(SI(1e-12,"s")) would compute the time development of the system for 12 pico seconds. For this example, we use the advance_time_ method simply to populate all the fields based on the initial values set for the magnetisation.

Extracting and saving data

We have three different ways of extracting data from the simulation

1. saving averaged values of fields (which can be analysed later)
2. saving spatiallialy resolved fields (which can be analysed later)
3. extracting field values at arbitrary positions from within the running program

In the sphere1 example discussed here, we use all three methods and will discuss these in more detail now:

Saving averaged data

The command sim.save_data_table() writes averages of all fields (see `Fields and subfields`_) into a text file. This file is best analysed using the ncol_ tool but can also just be read with a text editor. The format follows OOMMF’s odt fileformat: every row corresponds to one snap shot of the system (see save_data_table_).

The first and second line in the file are headers that explain the entity and the units of the entity saved in the corresponding column.

The ncol_ tool allows to extract particular columns easily so that these can be plotted later (useful for hysterises loop studies). In this example where have only one “timestep”, there is only one row of data in this file and we shall not explore this further here.

Extracting arbitrary data from the running programe

The line H_demag = sim.get_H_demag( [x,0,0] ) obtains the demagnetisation field at position (x,0,0). By default, the position is specified in SI units, and the data returned is also expressed in SI units.

The for-loop in the program (which changes x to range from -10*1e-9 to 10*1e-9 in steps of 1e-9) produces the following output

x = -1e-08 : H_demag =  None
x = -9e-09 : H_demag =  [-329655.76203912671, 130.62999726469423, 194.84338557811344]
x = -8e-09 : H_demag =  [-329781.46587966662, 66.963624669268853, 137.47161381890737]
x = -7e-09 : H_demag =  [-329838.57852402801, 181.46249265908259, 160.61298054099865]
x = -6e-09 : H_demag =  [-329899.63327447395, 131.06488858715838, 71.383139326493094]
x = -5e-09 : H_demag =  [-329967.79622912291, 82.209856975234786, -16.893046828024836]
x = -4e-09 : H_demag =  [-329994.67306536058, 61.622521557150371, -34.433041910642359]
x = -3e-09 : H_demag =  [-329997.62759666931, 23.222244635691535, -65.991127111463769]
x = -2e-09 : H_demag =  [-330013.90370482224, 10.11035370824321, -61.358763616681067]
x = -1e-09 : H_demag =  [-330023.50844056415, -6.9714476825652287, -54.900260456937708]
x = 0.0   : H_demag =  [-330030.98847923806, -26.808832466764223, -48.465748009067141]
x = 1e-09 : H_demag =  [-330062.38479507214, -38.660812022013424, -42.83439139610747]
x = 2e-09 : H_demag =  [-330093.78111090627, -50.512791577262625, -37.2030347831478]
x = 3e-09 : H_demag =  [-330150.72580001026, -64.552170478617398, -23.120555702674721]
x = 4e-09 : H_demag =  [-330226.19050178828, -77.236085707456397, -5.5373829923226916]
x = 5e-09 : H_demag =  [-330304.59300913941, -90.584413821813229, 14.090609104026118]
x = 6e-09 : H_demag =  [-330380.1392610991, -115.83746059068679, 37.072085708324757]
x = 7e-09 : H_demag =  [-330418.85831447819, -122.47512022500726, 62.379121138009992]
x = 8e-09 : H_demag =  [-330476.40747455234, -110.84257225592108, 108.06217226524763]
x = 9e-09 : H_demag =  [-330500.20126762061, -68.175725285038382, 162.46166752217249]
x = 1e-08 : H_demag =  [-330517.86675206106, -24.351273685146875, 214.40344001233677]

At position -1e-8, there is no field defined (this is just outside the mesh) and therefore the value None is returned.

We can see how the demagnetisation field varies slightly throughout the sphere. The x-component is approximately a third of the magnetisation, and the y- and z-components are close to zero (as would be expected for a perfectly round sphere).

Saving spatially resolved data (the fields)

The command sim.save_fields(all=True) will save all fields (see `Fields and subfields`_) spatially resolved for the current configuration into a file with name sphere1_dat.h5. The call sim.save_fields() will only save the magnetisation field (to store disk space). The data in this file is a compressed binary format (build on the hdf5_ standard) and can be extract and converted later using the nmagpp_ tool.

For example, to create a vtk_ file (for visualisation purposes) from the saved data, we can use:

nmagpp --vtk sphere1.vtk sphere1

where sphere1.vtk is the name of the vtk file that is to be generated.

Once this is executed, we can visualise the data. For this manual we use MayaVi_ as the visualisation tool for vtk files but there are others avialable (see vtk_).

We start MayaVi and load the vtk data file with mayavi -d sphere1.vtk. Using MayaVi’s menus, we can add a “VelocityVector” module to display the magnetisation:

The magnetisation is pointing in the x-direction (because we have set the initial magnetisation like this using the sim.set_m([1,0,0]) command).

The demagnetisation field should point in the opposite direction of the magnetisation. Let’s first plot a colour-coded plot of the value of the scalar magnetic potential phi from which the demag field is computed by taking the negative gradient:

We can see that the potential varies along the x-direction. The legend at the bottom of the figure shows the colour code used. We can also see from the title that the units of the potential phi are Ampere (this is the <A>).

Unless the user specifies a particular request for the units of data, the following rules apply:

• position are given in the same coordinates as the mesh coordinates (that is why the x, y and z axis have values going from -10 to 10).
• all field data is given in SI units.

The next plot shows the demag field (the vectors) together with isosurfaces of the magnetic potential:

It can be seen that the isosurfaces are completely flat planes (i.e. the potential is changing only along x) and the demagnetisation field is perpendicular to the isosurfaces. The colorbar on the left refers to the demagnetisation field which is expressed in Ampere per metre as can be seen from the label <A/m>.