more doc changes #252

doc/content/MSCA_EU_Horizon2020_Results/horizon2020_ex1.md

@@ -2,16 +2,13 @@
 
 # H2020-MSCA-IF-2019, Example: ground states
 
-!alert construction title=Documentation in-progress
-This section will have information that will be available when the review process is complete for the main publication encompassing this work. Please contact the developers if you need assistance with this aspect of the module.
-
 This page details how to obtain homogeneous solutions of BFO for the structural and spin order. The first calculation (`BFO_homogeneous_PA.i`) determines the polarization $\mathbf{P}$ and antiphase octahedral cage tilt vectors $\mathbf{A}$ along with the coupled elastic strain tensor components $\varepsilon_{ij}$ computed from the elastic displacement variable $\mathbf{u}$.
 
-The second (`BFO_homogeneous_LM.i`) uses information from the first to determine the antiferromagnetic order corresponding to the Neel vector $\mathbf{L}$ and total magnetization $\mathbf{m}$. We can separate these calculations because the magnetic system does not appreciably influence the relaxation dynamics of the polarization since the work required to move ions (i.e. those relative displacements corresponding to $\mathbf{P}$ and $\mathbf{A}$) is much larger than those to reorient the Fe spins.
+The second (`BFO_P0A0_mRD.i`) uses information from the first to determine the antiferromagnetic order corresponding to the Neel vector $\mathbf{L}$ and total magnetization $\mathbf{m}$. We can separate these calculations because the magnetic system does not appreciably influence the relaxation dynamics of the polarization since the work required to move ions (i.e. those relative displacements corresponding to $\mathbf{P}$ and $\mathbf{A}$) is much larger than those to reorient the Fe spins.
 
 # Homogeneous structural order
 
-We first consider a simple $2\times 2\times 2$ nm on a grid of $5 \times 5 \times 5$ finite elements with the `Mesh` block,
+We first consider a simple $1\times 1\times 1$ nm on a grid of $3 \times 3 \times 3$ finite elements with the `Mesh` block,
 
 !listing tutorial/BFO_homogeneous_PA.i
          block=Mesh
@@ -41,7 +38,7 @@ In this problem, the local strain is zero (no inhomogeneities) which means that
   \end{aligned}
 \end{equation}
 
-with $\Omega$ the computational volume. We find a set of global displacement vectors `disp_x, disp_y, disp_z` (see `AuxKernels`) such that the above condition is satisfied (see [!cite](Biswas2020) for an extended description of the method). In the input file, we see a number of objects that allow for this additional system,
+with $\Omega$ the computational volume. We find a set of global displacement vectors `disp_x, disp_y, disp_z` (see `AuxKernels`) such that the above condition is satisfied (see [!cite](Biswas2020) for an extended description of the method). In the input file, we see a number of MOOSE objects that allow for this additional system,
 
 !listing tutorial/BFO_homogeneous_PA.i
          block=ScalarKernels
@@ -79,10 +76,14 @@ where $Q_{ijkl}$ and $R_{ijkl}$ are the electrostrictive and rotostrictive coeff
          link=False
          language=python
 
-sets up the necessary residuals and jacobians for these equations. The initial conditions are set such that $\mathbf{P}$ and $\mathbf{A}$ are close to the expected values.
+sets up the necessary residuals and jacobians for these equations. The initial conditions are set such that $\mathbf{P}$ and $\mathbf{A}$ are close to the expected values. The boundary conditions (periodic) are both chosen to produce a homogeneous solution with $\mathbf{P}$ and $\mathbf{A}$ parallel ($\mathbf{P}\uparrow\uparrow\mathbf{A}$). We set this condition with the `BCs` block,
 
+!listing tutorial/BFO_homogeneous_PA.i
+         block=BCs
+         link=False
+         language=python
 
- and the boundary conditions (periodic) are both chosen to produce a homogeneous solution with $\mathbf{P}$ and $\mathbf{A}$ parallel ($\mathbf{P}\uparrow\uparrow\mathbf{A}$). The magnitudes of the structural order parameters are $P_s = |\mathbf{P}| = 0.945 \,\,\mathrm{C}/\mathrm{m}^2$  and $A_s = |\mathbf{A}| = 13.398^\circ$. The spontaneous $homogeneous$ normal and shear strain values are listed as rows with $\varepsilon_n = 1.308\times10^{-2}$ and $\varepsilon_s = 2.95 \times 10^{-3}$ respectively. The eight-fold domain variant symmetry is found in the pseudocubic reference frame corresponding to the below table.
+Note that the `DirichletBC` condition $\mathbf{u} = 0$ is set so there are no translations as $\mathbf{P}$ and $\mathbf{A}$ evolve to their minimum energy values. The resulting magnitudes of the structural order parameters are $P_s = |\mathbf{P}| = 0.945 \,\,\mathrm{C}/\mathrm{m}^2$  and $A_s = |\mathbf{A}| = 13.398^\circ$. The spontaneous $homogeneous$ normal and shear strain values are listed as rows with $\varepsilon_n = 1.308\times10^{-2}$ and $\varepsilon_s = 2.95 \times 10^{-3}$ respectively. The eight-fold domain variant symmetry is found in the pseudocubic reference frame corresponding to the below table.
 
 !equation
 \begin{array}{c|cccccccc}
@@ -98,15 +99,111 @@ sets up the necessary residuals and jacobians for these equations. The initial c
  \varepsilon_{xz}: & \varepsilon_s & -\varepsilon_s & -\varepsilon_s & \varepsilon_s & \varepsilon_s & -\varepsilon_s & -\varepsilon_s & \varepsilon_s \\
 \end{array}
 
-Note that $(\mathbf{P} \uparrow \downarrow \mathbf{A})$ is also a possible minimized energy solution of the thermodynamic potential $f$. Due to the symmetry of the electrostrictive and rotostrictive coupling terms, the table is left invariant under full reversal of $\mathbf{A}$. The free energy density of the eight-fold domain possibilities is -15.5653 $\mathrm{eV}\cdot\mathrm{nm}^{-3}$.
+Note that $(\mathbf{P} \uparrow \downarrow \mathbf{A})$ are also possible minimized energy solutions of the thermodynamic potential $f$. Due to the symmetry of the electrostrictive and rotostrictive coupling terms, the table is left invariant under full reversal of $\mathbf{A}$. The free energy density of the eight-fold domain possibilities is -15.5653 $\mathrm{eV}\cdot\mathrm{nm}^{-3}$.
 
+This simulation runs in 18.23 seconds on 4 processors.
 
 # Homogeneous spin order
 
-For this example, please consult `BFO_P0A0_mRD.i` in the tutorials subdirectory. We consider the output of the previous simulation $\{\mathbf{P},\mathbf{A}\}$ as an initial condition to find the homogeneous spin ground state. To do this, we use the following `Mesh` and `Variables` blocks, 
+For this example, please consult `BFO_P0A0_mRD.i` in the tutorials subdirectory. We consider the output of the previous simulation $\{\mathbf{P},\mathbf{A}\}$ as an initial condition to find the homogeneous spin ground state. To do this, we use the following `Mesh` block,
+
+!listing tutorial/BFO_P0A0_mRD.i
+         block=Mesh
+         link=False
+         language=python
+
+and `AuxVariables` blocks,
+
+!listing tutorial/BFO_P0A0_mRD.i
+         block=AuxVariables
+         link=False
+         language=python
+
+where we use (i.e. for vector component $P_z$) `initial_from_file_var = polar_z` and `initial_from_file_timestep = 'LATEST'` as code flags to load in the values of `polar_z` and others at the latest timestep from the Exodus output `BFO_P0A0.e`. In this sense, we have changed $\{\mathbf{P},\mathbf{A}\}$ to be an AuxVariable class instead of a Variable since they will be fixed for the duration of the evolution of the spins. Our variables to solve for are the set of sublattice magnetizations $\{\mathbf{m}_1, \mathbf{m}_2\}$ with postprocessed outputs $\mathbf{L} = \mathbf{m}_1 - \mathbf{m}_2$ and $\mathbf{m} = \mathbf{m}_1 + \mathbf{m}_2$. Note that the conventional factor of 2 is missing from the output of these calculations.
+
+This problem considers the evolution of the Landau-Lifshitz-Gilbert (Bloch) equation for the two sublattices,
 
-Our variables to solve for are the set of sublattice magnetizations $\{\mathbf{m}_1, \mathbf{m}_2\}$ with postprocessed outputs $\mathbf{L} = \mathbf{m}_1 - \mathbf{m}_2$ and $\mathbf{m} = \mathbf{m}_1 + \mathbf{m}_2$. Note that the conventional factor of 2 is missing from the output of these calculations.
+\begin{equation}\label{eqn:LLG_LLB}
+  \begin{aligned}
+    \frac{d \mathbf{m}_\eta}{dt} = -\frac{\gamma}{1+\alpha^2} \left(\mathbf{m}_\eta \times \mathbf{H}_\eta\right) - \frac{\gamma\alpha}{1+\alpha^2}\mathbf{m}_\eta\times\left(\mathbf{m}_\eta\times\mathbf{H}_\eta\right) + \frac{\gamma\tilde\alpha_\parallel}{(1+\alpha^2)} m_{\eta}^2 \left[ m_{\eta}^2 - 1 \right]\mathbf{m}_\eta,
+  \end{aligned}
+\end{equation}
+
+where $\gamma$ corresponds to the electron gyromagnetic factor equal to $2.2101\times10^5$ rad. m A${}^{-1}$ s ${}^{-1}$, $\alpha$ the Gilbert damping, and $\tilde\alpha_\parallel$ the longitudinal damping constant from the LLB approximation. The effective fields are defined as $\mathbf{H}_\eta = - \mu_0^{-1} M_s^{-1} \delta f_\mathrm{mag} / \delta \mathbf{m}_\eta$ with $\mu_0$ the permeability of vacuum. The saturation magnetization density of the BFO sublattices is $M_s = 4.0$ $\mu$B/Fe [!cite](Dixit2015). We choose an initial condition close to the expected ground state values. We set $\alpha = 0.1$ and $\tilde\alpha_\parallel = 10^3$ and ringdown the magnetic system to an energy minimum.
+
+The spin free energy, $f_\mathrm{mag}$, is correspondingly,
+
+\begin{equation}\label{eqn:fmag}
+  \begin{aligned}
+    f_\mathrm{mag} = f_\mathrm{exch} + f_\mathrm{DMI} + f_\mathrm{easy} + f_\mathrm{anis}
+  \end{aligned}
+\end{equation}
+
+where
+
+\begin{equation}\label{eqn:fmag_split}
+  \begin{aligned}
+    f_\mathrm{exch} &= D_e \left(\mathbf{L}^2 - \mathbf{m}^2\right)\\
+    f_\mathrm{DMI} &= 4 D_0 \mathbf{A} \cdot \left( \mathbf{L} \times \mathbf{m}\right) \\
+    f_\mathrm{easy} &= \sum\limits_{\eta = 1}^2 K_1 \left(\textbf{m}_\eta \cdot \hat{\mathbf{P}}\right)^2 \\
+    f_\mathrm{anis} &= \sum\limits_{\eta = 1}^2 \left(K_1^c + a|\mathbf{A}|^2\right)\left(m_{\eta,x}^2 m_{\eta,y}^2 m_{\eta,z}^2\right).
+  \end{aligned}
+\end{equation}
+
+Each of these contributions to the effective fields $\mathbf{H}$ contain different physics. For example, the effective fields from $f_\mathrm{exch}$ competes with those from $f_\mathrm{DMI}$ to provide a nearly-collinear spin structure with a weak canting angle $\phi^\mathrm{WFM} = \cos^{-1}{\left(\mathbf{m}_1\cdot \mathbf{m}_2\right)}$. The introduction of an effective field due to $f_\mathrm{easy}$ puts the collinearity within the easy plane defined by the director $\hat{\mathbf{P}}$. This yields $\theta_1 = \theta_2 = 90^\circ$ where $\theta_\eta = \cos^{-1}{\left(\mathbf{m}_\eta \cdot \hat{\mathbf{P}}\right)}$. The weak effective field generated from $f_\mathrm{anis}$ splits the easy-plane degeneracy into a six-fold symmetry upon energy minimization. The `Materials` block,
+
+!listing tutorial/BFO_P0A0_mRD.i
+         block=Materials
+         link=False
+         language=python
+
+seeds the values of these quantities for $f_\mathrm{mag}$ (note that they are set at the top of the input file for ease of use). This is because the MOOSE `Parser` reads in `${quantity}` from the out-of-block definition `quantity = something`. The definition for `alpha_long` corresponds to that of $\tilde\alpha_\parallel$ and we set it with a constant `ParsedFunction`. It could be defined the same way as the other materials coefficients but we leave this as an option for the user in case they want to give $\tilde\alpha_\parallel$ some functional dependence (on time or space for example).
+
+The `Kernels` block is used to initialize the partial differential equation and computes the residual and jacobian contributions for each component of the sublattices,
+
+!listing tutorial/BFO_P0A0_mRD.i
+         block=Kernels
+         link=False
+         language=python
+
+One can see the derivations of the residual and jacobian entries for the `Kernels` `Objects` by visiting the `Syntax` page. Some more options for this problem are located in the `Executioner` block (i.e. time step size `dt`),
+
+!listing tutorial/BFO_P0A0_mRD.i
+         block=Executioner
+         link=False
+         language=python
+
+where we have used the [NewmarkBeta](https://mooseframework.inl.gov/source/timeintegrators/NewmarkBeta.html) time integration method. We find that this is most numerically stable for AFM ringdown problems. The simulation runs fairly quickly (400 seconds) on 4 processors, stepping through a few thousands of time steps to ringdown $\mathbf{m}_\eta$. A visualization of the components of $\mathbf{m}_\eta$ from the `ParaView` filter `PlotGlobalVariablesOverTime` is provided in the below figure.
+
+!media media/tut_AFM_BFO.png style=display:block;margin:auto;width:60%; caption=Angular quantities $\phi^\mathrm{WFM}, \theta_1,$ and $\theta_2$ during ringdown. id=tut_AFM_ground
+
+By setting the initial condition of $\mathbf{m}_\eta$ in different directions, one can obtain the below table for the eight-fold orientations of $\mathbf{P}\uparrow \mathbf{A}$.
+
+!equation
+\begin{array}{c|cccccccc}
+\hline
+\hline
+ & & & & & & & & \\
+ (\mathbf{P}\uparrow\uparrow\mathbf{A}): & [111] & [\bar{1}11] & [1\bar{1}\bar{1}] & [\bar{1}\bar{1}\bar{1}] & [1\bar{1}1] & [11\bar{1}] & [\bar{1}\bar{1}1] & [\bar{1}1\bar{1}] \\
+ \hline
+ & & & & & & & & \\
+\mathbf{L} \simeq & [\bar{1}10] & [101] & [101] & [\bar{1}10] & [110] & [\bar{1}10] & [\bar{1}10] & [110] \\
+& [\bar{1}01] & [110] & [110] & [\bar{1}01] & [011] & [011] & [011] & [011] \\
+& [0\bar{1}1] & [01\bar{1}] & [01\bar{1}] & [0\bar{1}1] & [\bar{1}01] & [101] & [101] & [\bar{1}01] \\
+& [1\bar{1}0] & [\bar{1}0\bar{1}] & [\bar{1}0\bar{1}] & [1\bar{1}0]  & [\bar{1}\bar{1}0] & [1\bar{1}0] & [1\bar{1}0] & [0\bar{1}\bar{1}] \\
+& [10\bar{1}] & [\bar{1}\bar{1}0] & [\bar{1}\bar{1}0] & [10\bar{1}] & [0\bar{1}\bar{1}] & [0\bar{1}\bar{1}] & [0\bar{1}\bar{1}] & [0\bar{1}\bar{1}] \\
+& [01\bar{1}] & [0\bar{1}1] & [0\bar{1}1] & [01\bar{1}] & [10\bar{1}] & [\bar{1}0\bar{1}] & [\bar{1}0\bar{1}] & [10\bar{1}] \\
+\hline
+& & & & & & & & & & & & & & & & \\
+\mathbf{m} \simeq & [\bar{1}\bar{1}2] & [12\bar{1}] & [\bar{1}\bar{2}1] & [11\bar{2}] & [\bar{1}12] & [112] & [\bar{1}\bar{1}\bar{2}] & [1\bar{1}\bar{2}] \\
+ & [1\bar{2}1] & [\bar{1}1\bar{2}] & [1\bar{1}2] & [\bar{1}2\bar{1}] & [\bar{2}\bar{1}1] & [2\bar{1}1] & [\bar{2}1\bar{1}] & [21\bar{1}] \\
+ & [2\bar{1}\bar{1}] & [\bar{2}\bar{1}\bar{1}] & [211] & [\bar{2}11] & [\bar{1}\bar{2}\bar{1}] & [1\bar{2}\bar{1}] & [\bar{1}21] & [121]\\
+ & [11\bar{2}] & [\bar{1}\bar{2}1] & [12\bar{1}] & [\bar{1}\bar{1}2] & [1\bar{1}\bar{2}] & [\bar{1}\bar{1}\bar{2}] & [112] & [\bar{1}12] \\
+ & [\bar{1}2\bar{1}] & [1\bar{1}2] & [\bar{1}1\bar{2}] & [1\bar{2}1] & [21\bar{1}] & [\bar{2}1\bar{1}] & [2\bar{1}1] & [\bar{2}\bar{1}1] \\
+ & [\bar{2}11] & [211] & [\bar{2}\bar{1}\bar{1}] & [2\bar{1}\bar{1}] & [121] & [\bar{1}21] & [1\bar{2}\bar{1}] & [\bar{1}\bar{2}\bar{1}] \\
+\end{array}
 
+Importantly, the nonzero value of $\mathbf{m}$ corresponds to a canted angle $\phi^\mathrm{WFM}$ of about $1.21^\circ$ which agrees well with the literature on BFO.
 
 
 This project [SCALES - 897614](https://cordis.europa.eu/project/id/897614) was funded for 2021-2023 at the [Luxembourg Institute of Science and Technology](https://www.list.lu/) under principle investigator [Jorge Íñiguez](https://sites.google.com/site/jorgeiniguezresearch/). The research was carried out within the framework of the [Marie Skłodowska-Curie Action (H2020-MSCA-IF-2019)](https://ec.europa.eu/info/funding-tenders/opportunities/portal/screen/opportunities/topic-details/msca-if-2020) fellowship.

doc/content/MSCA_EU_Horizon2020_Results/horizon2020_results2.md

@@ -44,7 +44,7 @@ where $Q_{ijkl}$ and $R_{ijkl}$ are the electrostrictive and rotostrictive coeff
 
 Note that $(\mathbf{P} \uparrow \downarrow \mathbf{A})$ is also a possible minimized energy solution of the thermodynamic potential $f$. Due to the symmetry of the electrostrictive and rotostrictive coupling terms, the table is left invariant under full reversal of $\mathbf{A}$. The free energy density of the eight-fold domain possibilities is -15.5653 $\mathrm{eV}\cdot\mathrm{nm}^{-3}$.
 
-To reproduce these results, we refer the reader to our examples documentation [here](MSCA_EU_Horizon2020_Results/horizon2020_ex1.md).
+To reproduce these results, we refer the reader to our examples documentation [here](MSCA_EU_Horizon2020_Results/horizon2020_ex1.md) for more details on the input file.
 
 # Antiferromagnetic ringdown (homogeneous spin state)
 
@@ -90,7 +90,7 @@ By selecting initial conditions corresponding to the six-fold possible orientati
  & [\bar{2}11] & [211] & [\bar{2}\bar{1}\bar{1}] & [2\bar{1}\bar{1}] & [121] & [\bar{1}21] & [1\bar{2}\bar{1}] & [\bar{1}\bar{2}\bar{1}] \\
 \end{array}
 
-To reproduce these results, we refer the reader to our examples documentation [here](MSCA_EU_Horizon2020_Results/horizon2020_ex1.md).
+To reproduce these results, we refer the reader to our examples documentation [here](MSCA_EU_Horizon2020_Results/horizon2020_ex1.md) for details on the input file.
 
 # Structural domain walls
 

tutorial/BFO_P0A0_mRD.i

@@ -1,10 +1,11 @@
 
-Dedef = 3.0
-D0def = 0.01
-K1def = -0.2
-K1cdef = -0.01
-Ktdef = -0.0001
-alphadef = 0.01
+Dedef = 3.7551
+D0def = 0.003
+K1def = -5.0068
+K1cdef = -0.0550748
+Ktdef = -0.00365997
+
+alphadef = 0.1
 
 [Mesh]
   [fileload]
@@ -584,6 +585,22 @@ alphadef = 0.01
     type = TimestepSize
   [../]
 
+  [./ph]
+    type = ElementAverageValue
+    execute_on = 'initial timestep_end final'
+    variable = ph
+  [../]
+  [./th1]
+    type = ElementAverageValue
+    execute_on = 'initial timestep_end final'
+    variable = th1
+  [../]
+  [./th2]
+    type = ElementAverageValue
+    execute_on = 'initial timestep_end final'
+    variable = th2
+  [../]
+
   [./FafmSLexch]
     type = AFMSublatticeSuperexchangeEnergy
     execute_on = 'initial timestep_end final'
@@ -671,12 +688,19 @@ alphadef = 0.01
     dt = dt 
     execute_on = 'timestep_end final'
   [../]
+
+
+  [./elapsed]
+    type = PerfGraphData
+    section_name = "Root"  # for profiling the problem
+    data_type = total
+  [../]
 []
 
 [UserObjects]
   [./kill]
     type = Terminator
-    expression = 'perc_change <= 1.0e-8'
+    expression = 'perc_change <= 1.0e-5'
   [../]
 []
 
@@ -698,16 +722,14 @@ alphadef = 0.01
   [../]
 
   dtmin = 1e-18
-  dtmax = 1.25e-7
+  dtmax = 1.0e-6
 
   [./TimeStepper]
     type = IterationAdaptiveDT
     optimal_iterations = 10  #usually 10
     linear_iteration_ratio = 100
     dt = 1e-8
   [../]
-
-  num_steps = 20
 []
 
 [Outputs]
@@ -716,7 +738,7 @@ alphadef = 0.01
     type = Exodus
     file_base = out_BFO_P0A0_mRD
     elemental_as_nodal = true
-    interval = 2
+    interval = 4
   [../]
 []