#!/usr/bin/env python
# coding: utf-8

# # Micromagnetic standard problem 4
# 
# ## Problem specification
# 
# The sample is a thin film cuboid with dimensions:
# 
# - length $l_{x} = 500 \,\text{nm}$,
# - width $l_{y} = 125 \,\text{nm}$, and
# - thickness $l_{z} = 3 \,\text{nm}$.
# 
# The material parameters (similar to permalloy) are:
# 
# - exchange energy constant $A = 1.3 \times 10^{-11} \,\text{J/m}$,
# - magnetisation saturation $M_\text{s} = 8 \times 10^{5} \,\text{A/m}$.
# 
# Magnetisation dynamics are governed by the Landau-Lifshitz-Gilbert equation
# 
# $$\frac{d\mathbf{m}}{dt} = \underbrace{-\gamma_{0}(\mathbf{m} \times \mathbf{H}_\text{eff})}_\text{precession} + \underbrace{\alpha\left(\mathbf{m} \times \frac{d\mathbf{m}}{dt}\right)}_\text{damping}$$
# 
# where $\gamma_{0} = 2.211 \times 10^{5} \,\text{m}\,\text{A}^{-1}\,\text{s}^{-1}$ and Gilbert damping $\alpha=0.02$.
# 
# In the standard problem 4, the system is first relaxed at zero external magnetic field and then, starting from the obtained equlibrium configuration, the magnetisation dynamics are simulated for two external magnetic fields $\mathbf{B}_{1} = (-24.6, 4.3, 0.0) \,\text{mT}$ and $\mathbf{B}_{2} = (-35.5, -6.3, 0.0) \,\text{mT}$.
# 
# More detailed specification of Standard problem 4 can be found in Ref. 1.
# 
# ## Simulation
# 
# 
# In the first step, we import the required `discretisedfield` and `oommfc` modules.

# In[1]:


import discretisedfield as df
import oommfc as oc


# Now, we can set all required geometry and material parameters.

# In[2]:


# Geometry
lx = 500e-9  # x dimension of the sample(m)
ly = 125e-9  # y dimension of the sample (m)
lz = 3e-9  # sample thickness (m)

# Material (permalloy) parameters
Ms = 8e5  # saturation magnetisation (A/m)
A = 1.3e-11  # exchange energy constant (J/m)

# Dynamics (LLG equation) parameters
gamma = 2.211e5  # gyromagnetic ratio (m/As)
alpha = 0.02  # Gilbert damping


# ## First stage
# 
# In the first stage, we need to relax the system at zero external magnetic field.
# 
# We choose `stdprob4` to be the name of the system. This name will be used to name all output files created by OOMMF.

# In[3]:


system = oc.System(name='stdprob4')


# In order to completely define the micromagnetic system, we need to provide:
# 
# 1. hamiltonian $\mathcal{H}$
# 2. dynamics $\text{d}\mathbf{m}/\text{d}t$
# 3. magnetisation $\mathbf{m}$
# 
# The mesh is created by providing two points `p1` and `p2` between which the mesh domain spans and the size of a discretisation cell. We choose the discretisation to be $(5, 5, 3) \,\text{nm}$.

# In[4]:


cell = (5e-9, 5e-9, 3e-9)  # mesh discretisation (m)

mesh = oc.Mesh(p1=(0, 0, 0), p2=(lx, ly, lz), cell=cell)  # Create a mesh object.


# We can visualise the mesh domain and a discretisation cell:

# In[5]:


get_ipython().run_line_magic('matplotlib', 'inline')
mesh


# **Hamiltonian:** In the second step, we define the system's Hamiltonian. In this standard problem, the Hamiltonian contains only exchange and demagnetisation energy terms. Please note that in the first simulation stage, there is no applied external magnetic field. Therefore, we do not add Zeeman energy term to the Hamiltonian.

# In[6]:


system.hamiltonian = oc.Exchange(A) + oc.Demag()


# We can check what is the continuous model of system's Hamiltonian.

# In[7]:


system.hamiltonian


# **Dynamics:** Similarly, the system's dynamics is defined by providing precession and damping terms (LLG equation).

# In[8]:


system.dynamics = oc.Precession(gamma) + oc.Damping(alpha)

system.dynamics  # check the dynamics equation


# **Magnetisation:** Finally, we have to provide the magnetisation configuration that is going to be relaxed subsequently. We choose the uniform configuration in $(1, 0.25, 0.1)$ direction, and as norm (magnitude) we set the magnetisation saturation $M_\text{s}$. In order to create the magnetisation configuration, we create a `Field` object from the `discretisedfield` module.

# In[9]:


system.m = df.Field(mesh, value=(1, 0.25, 0.1), norm=Ms)


# Now, the system is fully defined.
# 
# **Energy minimisation:** The system (its magnetisation) is evolved using a particular driver. In the first stage, we need to relax the system - minimise its energy. Therefore, we create `MinDriver` object and drive the system using its `drive` method.

# In[10]:


md = oc.MinDriver()  # create energy minimisation driver
md.drive(system)  # minimise the system's energy


# The system is now relaxed. We can now obtain some data characteristic to the magnetisation field.

# In[11]:


print('The average magnetisation is {}.'.format(system.m.average))

print('The magnetisation at the mesh centre {} is {}.'.format(
        system.m.mesh.centre, system.m(system.m.mesh.centre)))


# ## Second stage: field $\mathbf{B}_{1}$

# In the second stage, we need to apply an external magnetic field $\mathbf{B}_{1} = (-24.6, 4.3, 0.0) \,\text{mT}$ to the system. In other words, we have to add Zeeman energy term to the Hamiltonian.

# In[12]:


# Add Zeeman energy term to the Hamiltonian
H1 = (-24.6e-3/oc.mu0, 4.3e-3/oc.mu0, 0.0)
system.hamiltonian += oc.Zeeman(H1)


# If we now inspect the Hamiltonian, we see that an additional Zeeman term is added.

# In[13]:


system.hamiltonian


# Finally, we can run the simulation using `TimeDriver` this time. We run the magnetisation evolution for $t=1 \,\text{ns}$, during which we save the system's state $n=200$ times.

# In[14]:


t = 1e-9  # simulation time (s)
n = 200  # number of data saving steps

td = oc.TimeDriver()  # create time driver
td.drive(system, t=t, n=n)  # drive the system


# ### Postprocessing

# When we drove the system using the `TimeDriver`, we specified that we want to save the magnetisation configuration $n=200$ times. A detailed table of all computed parameters from the last simulation can be shown from the datatable (`system.dt`), which is a `pandas` dataframe [2].
# 
# For instance, if we want to show the last 10 rows in the table, we run:

# In[15]:


system.dt.tail()


# Finally, we want to plot the average magnetisation configuration `my` as a function of time `t`:

# In[16]:


myplot = system.dt.plot("t", "my")


# ## References
# 
# [1] µMAG Site Directory: http://www.ctcms.nist.gov/~rdm/mumag.org.html
#     
# [2] Pandas: http://pandas.pydata.org/