diff --git a/.github/workflows/docker-image.yml b/.github/workflows/docker-image.yml
index e9d1335..a002cd2 100644
--- a/.github/workflows/docker-image.yml
+++ b/.github/workflows/docker-image.yml
@@ -17,6 +17,8 @@ jobs:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v3
+ with:
+ submodules: 'true'
- name: Login to GH container registry
uses: docker/login-action@v1
with:
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 6ef1a31..dd9d2a9 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -30,6 +30,8 @@ jobs:
steps:
# Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
- uses: actions/checkout@v3
+ with:
+ submodules: 'true'
- name: Python wheels manylinux stable build
uses: RalfG/python-wheels-manylinux-build@v0.5.0
with:
@@ -53,7 +55,9 @@ jobs:
os: [macos-latest, windows-latest]
python-version: ['3.8', '3.9', '3.10', '3.11', '3.12']
steps:
- - uses: actions/checkout@v2
+ - uses: actions/checkout@v3
+ with:
+ submodules: 'true'
- name: Set up Python
uses: actions/setup-python@v2
with:
diff --git a/.github/workflows/python-test.yml b/.github/workflows/python-test.yml
index f4431b4..b784570 100644
--- a/.github/workflows/python-test.yml
+++ b/.github/workflows/python-test.yml
@@ -34,6 +34,8 @@ jobs:
steps:
- uses: actions/checkout@v3
+ with:
+ submodules: 'true'
- name: Set up Python ${{ matrix.python-version }}
uses: actions/setup-python@v3
with:
@@ -57,6 +59,8 @@ jobs:
steps:
- uses: actions/checkout@v3
+ with:
+ submodules: 'true'
- name: Set up Python ${{ matrix.python-version }}
uses: actions/setup-python@v3
with:
diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..5f164bb
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "third_party/kissfft"]
+ path = third_party/kissfft
+ url = https://github.com/mborgerding/kissfft
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 69aec4e..e3d95f8 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -31,6 +31,7 @@ repos:
- id: prettier
files: \.(html|json|markdown|md|yaml|yml)$
exclude: (^docs/api/.*)
+
- repo: https://github.com/pycqa/isort
rev: 5.10.1
hooks:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6136684..ff4115d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,14 @@
# Changelog
+# 1.5.0 (WIP)
+
+- Dipole interaction added to the `SB Model`
+- Kasdin 1/f noise generator added to the `noise` module and to the solvers
+- reworking the solvers for better performance and stability
+- added a simple noise model to the `utils` class. It exists outside standard simulation procedures.
+- added LLGB bindings and code. The solver is still WIP and doesn't integrate with more advanced features yet.
+- added aliases for `ScalarDriver` -- for example, instead of calling `ScalarDriver.getConstantDriver`, you can now call `constantDriver` directly to create a constant driver.
+
# 1.4.1
- Adding a basic optimisation script in the `optimization` module.
diff --git a/README.md b/README.md
index 64e700c..e3d4b29 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-
+
# CMTJ
@@ -9,15 +9,34 @@
[![Streamlit](https://static.streamlit.io/badges/streamlit_badge_black_white.svg)](http://cmtj-simulations.streamlit.app/)
![Downloads](https://img.shields.io/pypi/dm/cmtj.svg)
+## Table of contents
+
+- [CMTJ](#cmtj)
+ - [Table of contents](#table-of-contents)
+ - [Short description](#short-description)
+ - [Web GUI](#web-gui)
+ - [Quickstart](#quickstart)
+ - [Installation :rocket:](#installation-rocket)
+ - [Extra dependencies](#extra-dependencies)
+ - [Documentation and examples](#documentation-and-examples)
+ - [Extensions](#extensions)
+ - [Citing](#citing)
+- [Development](#development)
+ - [Acknowledgements](#acknowledgements)
+ - [Contributions](#contributions)
+ - [Docker](#docker)
+ - [Precommit](#precommit)
+ - [Documentation builds](#documentation-builds)
+
## Short description
The `cmtj` name may be misleading -- the MTJ (Magnetic Tunnel Junctions) are not the only structures that may be simulated.
The library allows for macromagnetic simulation of various multilayer spintronic structures. The package uses C++ implementation of (s)LLGS (stochastic Landau-Lifschitz-Gilbert-Slonczewski) equation with various field contributions included for instance: anisotropy, interlayer exchange coupling, demagnetisation, dipole fields etc.
It is also possible to connect devices in parallel or in series to have electrically coupled arrays.
-## Demo
+## Web GUI
-Check out the [streamlit hosted demo here](http://cmtj-simulations.streamlit.app/).
+Check out the [streamlit hosted demo here](http://cmtj-simulations.streamlit.app/). You can simulate PIMM spectra and Spin-Diode spectra there. Let us know if you have any issues with the demo.
## Quickstart
@@ -63,31 +82,10 @@ The package requires (if `utils` subpackage is used):
- matplotlib
```
-## Subpackages
-```mermaid
-graph TD;
-cmtj --> models
-cmtj --> utils
-models --> domain_dynamics
-models --> drivers
-models --> ensemble
-models --> general_sb
-models --> oersted
-utils --> optimization
-utils --> parallel
-utils --> plotting
-utils --> procedures
-utils --> resistance
-utils --> solvers
-utils --> linear
-utils --> energy
-utils --> filters
-```
-
-
-## Read the docs
+## Documentation and examples
Documentation: [https://lemurpwned.github.io/cmtj](https://lemurpwned.github.io/cmtj)
+There are many examples available, check out the [examples section in the docs](https://lemurpwned.github.io/cmtj/experimental-methods/introduction/)
## Extensions
@@ -155,10 +153,19 @@ There are couple of stages to building the documentation
This is mostly for the C++ documentation. Furture changes may couple C++ and Python docs.
2. Build stubs
The stubgen is `pybind11-stubgen` or `mypy stubgen` with the latter being preferred now.
- E.g. to generate `Stack` module stubs we can go:
+ Before running the stubgen, make sure to install the package with:
+ ```
+ python3 -m pip install .
+ ```
+ avoid using `-e` flag as it may cause issues with the stubgen.
+ Then to generate, for instance, `Stack` module stubs we can do:
```
stubgen -m cmtj.stack -o target-stub-dir/
```
+ or
+ ```
+ python3 -c "import mypy.stubgen; mypy.stubgen.main(['-p', 'cmtj.stack', '-o', 'target-stub-dir/'])"
+ ```
More info here: https://mypy.readthedocs.io/en/stable/stubgen.html.
3. Parse stubs to Markdown.
This stage is done by running:
diff --git a/cmtj/__init__.pyi b/cmtj/__init__.pyi
index c792583..640c347 100644
--- a/cmtj/__init__.pyi
+++ b/cmtj/__init__.pyi
@@ -1,136 +1,195 @@
+import typing
from typing import Any, ClassVar, Dict, List, overload
xaxis: Axis
yaxis: Axis
zaxis: Axis
-
+none: Axis
+all: Axis
def c_dot(arg0: CVector, arg1: CVector) -> float:
"""Compute dot (scalar) product of two CVectors."""
...
+def constantDriver(constant: float) -> ScalarDriver:
+ """
+ Constant driver produces a constant signal of a fixed amplitude.
+ :param constant: constant value of the driver (constant offset/amplitude)
+ """
+ ...
-class AxialDriver:
- @overload
- def __init__(self, x_driver: ScalarDriver, y_driver: ScalarDriver,
- z_driver: ScalarDriver) -> None:
- ...
+def sineDriver(
+ constantValue: float, amplitude: float, frequency: float, phase: float
+) -> ScalarDriver:
+ """
+ Produces a sinusoidal signal with some offset (constantValue), amplitude frequency and phase offset.
+ :param constantValue: vertical offset. The sine will oscillate around this value.
+ :param amplitude: amplitude of the sine wave
+ :param frequency: frequency of the sine
+ :param phase: phase of the sine in radians.
+ """
+ ...
- @overload
- def __init__(self, arg0: List[ScalarDriver]) -> None:
- ...
+def gaussianImpulseDriver(
+ constantValue: float, amplitude: float, t0: float, sigma: float
+) -> ScalarDriver:
+ """
+ Gaussian impulse driver. It has amplitude starts at t0 and falls off with sigma.
- @overload
- def __init__(*args, **kwargs) -> Any:
- ...
+ Formula:
+ A * exp(-((t - t0) ** 2) / (2 * sigma ** 2))
- @overload
- def applyMask(self, arg0: CVector) -> None:
- ...
+ :param constantValue: offset of the pulse (vertical)
+ :param amplitude: amplitude that is added on top of the constantValue
+ :param t0: start of the pulse
+ :param sigma: fall-off of the Gaussian pulse
+ """
+ ...
- @overload
- def applyMask(self, arg0: List[int]) -> None:
- ...
+def gaussianStepDriver(
+ constantValue: float, amplitude: float, t0: float, sigma: float
+) -> ScalarDriver:
+ """Gaussian step driver (erf function). It has amplitude starts at t0 and falls off with sigma.
- @overload
- def applyMask(*args, **kwargs) -> Any:
- ...
+ Formula:
+ f(t) = constantValue + amplitude * (1 + erf((t - t0) / (sigma * sqrt(2))))
- def getCurrentAxialDrivers(self, arg0: float) -> CVector:
- ...
+ :param constantValue: offset of the pulse (vertical)
+ :param amplitude: amplitude that is added on top of the constantValue
+ :param t0: start of the pulse
+ :param sigma: fall-off of the Gaussian pulse
+ """
+ ...
- def getVectorAxialDriver(self, arg0: float, arg1: float) -> AxialDriver:
- ...
+def posSineDriver(
+ constantValue: float, amplitude: float, frequency: float, phase: float
+) -> ScalarDriver:
+ """Produces a positive sinusoidal signal with some offset (constantValue), amplitude frequency and phase offset.
+ :param constantValue: vertical offset. The sine will oscillate around this value.
+ :param amplitude: amplitude of the sine wave
+ :param frequency: frequency of the sine
+ :param phase: phase of the sine in radians.
+ """
+ ...
+
+def pulseDriver(
+ constantValue: float, amplitude: float, period: float, cycle: float
+) -> ScalarDriver:
+ """
+ Produces a square pulse of certain period and cycle
+ :param constantValue: offset (vertical) of the pulse. The pulse amplitude will be added to this.
+ :param amplitude: amplitude of the pulse signal
+ :param period: period of the signal in seconds
+ :param cycle: duty cycle of the signal -- a fraction between [0 and 1].
+ """
+ ...
+
+def stepDriver(
+ constantValue: float, amplitude: float, timeStart: float, timeStop: float
+) -> ScalarDriver:
+ """
+ Get a step driver. It has amplitude between timeStart and timeStop and 0 elsewhere
+ :param constantValue: offset of the pulse (vertical)
+ :param amplitude: amplitude that is added on top of the constantValue
+ :param timeStart: start of the pulse
+ :param timeStop: when the pulse ends
+ """
+ ...
+def trapezoidDriver(
+ constantValue: float,
+ amplitude: float,
+ timeStart,
+ edgeTime: float,
+ steadyTime: float,
+) -> ScalarDriver:
+ """Create Trapezoid driver. Has a rising and a falling edge.
+ :param constantValue: offset of the pulse (vertical)
+ :param amplitude: amplitude that is added on top of the constantValue
+ :param timeStart: start of the pulse
+ :param edgeTime: time it takes to reach the maximum amplitude
+ :param steadyTime: time it spends in a steady state
+ """
+ ...
+
+class AxialDriver:
+ @overload
+ def __init__(
+ self, x_driver: ScalarDriver, y_driver: ScalarDriver, z_driver: ScalarDriver
+ ) -> None: ...
+ @overload
+ def __init__(self, arg0: List[ScalarDriver]) -> None: ...
+ @overload
+ def __init__(*args, **kwargs) -> Any: ...
+ @overload
+ def applyMask(self, arg0: CVector) -> None: ...
+ @overload
+ def applyMask(self, arg0: List[int]) -> None: ...
+ @overload
+ def applyMask(*args, **kwargs) -> Any: ...
+ def getCurrentAxialDrivers(self, arg0: float) -> CVector: ...
+ def getVectorAxialDriver(self, arg0: float, arg1: float) -> AxialDriver: ...
class Axis:
__entries: Any = ...
xaxis: Any = ...
yaxis: Any = ...
zaxis: Any = ...
-
- def __init__(self, value: int) -> None:
- ...
-
- def __eq__(self, other: object) -> bool:
- ...
-
- def __getstate__(self) -> int:
- ...
-
- def __hash__(self) -> int:
- ...
-
- def __index__(self) -> int:
- ...
-
- def __int__(self) -> int:
- ...
-
- def __ne__(self, other: object) -> bool:
- ...
-
- def __setstate__(self, state: int) -> None:
- ...
-
+ none: Any = ...
+ all: Any = ...
+
+ def __init__(self, value: int) -> None: ...
+ def __eq__(self, other: object) -> bool: ...
+ def __getstate__(self) -> int: ...
+ def __hash__(self) -> int: ...
+ def __index__(self) -> int: ...
+ def __int__(self) -> int: ...
+ def __ne__(self, other: object) -> bool: ...
+ def __setstate__(self, state: int) -> None: ...
@property
- def name(self) -> Any:
- ...
-
+ def name(self) -> Any: ...
@property
- def __doc__(self) -> Any:
- ...
-
+ def __doc__(self) -> Any: ...
@property
- def __members__(self) -> Any:
- ...
-
+ def __members__(self) -> Any: ...
class CVector:
- def __init__(self, x: float, y: float, z: float) -> None:
- ...
-
- def length(self) -> float:
- ...
-
+ def __init__(self, x: float, y: float, z: float) -> None: ...
+ def length(self) -> float: ...
+ def normalize(self) -> None: ...
+ def tolist(self) -> List[float]: ...
+ def __add__(self, arg0: CVector) -> CVector: ...
+ def __eq__(self, arg0: CVector) -> bool: ...
+ def __getitem__(self, arg0: int) -> float: ...
+ def __iter__(self) -> typing.Iterator[float]: ...
+ def __iadd__(self, arg0: CVector) -> CVector: ...
+ def __imul__(self, arg0: float) -> CVector: ...
+ def __isub__(self, arg0: CVector) -> CVector: ...
+ def __len__(self) -> int: ...
+ def __mul__(self, arg0: float) -> CVector: ...
+ def __ne__(self, arg0: CVector) -> bool: ...
+ def __rmul__(self, arg0: float) -> CVector: ...
+ def __sub__(self, arg0: CVector) -> CVector: ...
@property
- def x(self) -> float:
- ...
-
+ def x(self) -> float: ...
@x.setter
- def x(self, val: float) -> None:
- ...
-
+ def x(self, val: float) -> None: ...
@property
- def y(self) -> float:
- ...
-
+ def y(self) -> float: ...
@y.setter
- def y(self, val: float) -> None:
- ...
-
+ def y(self, val: float) -> None: ...
@property
- def z(self) -> float:
- ...
-
+ def z(self) -> float: ...
@z.setter
- def z(self, val: float) -> None:
- ...
-
+ def z(self, val: float) -> None: ...
class Junction:
@overload
- def __init__(self, layers: List[Layer], filename: str = ...) -> None:
- ...
-
+ def __init__(self, layers: List[Layer], filename: str = ...) -> None: ...
@overload
- def __init__(self,
- layers: List[Layer],
- filename: str,
- Rp: float = ...,
- Rap: float = ...) -> None:
- ...
-
+ def __init__(
+ self, layers: List[Layer], filename: str, Rp: float = ..., Rap: float = ...
+ ) -> None: ...
@overload
def __init__(
self,
@@ -149,37 +208,31 @@ class Junction:
length of the layers passed (they directly correspond to each layer).
Calculates the magnetoresistance as per: __see reference__:
Spin Hall magnetoresistance in metallic bilayers by Kim, J. et al.
- :param Rx0
- :param Ry0
- :param AMR_X
- :param AMR_Y
- :param SMR_X
- :param SMR_Y
- :param AHE
+ :param Rx0: Magnetoresistance offset longitudinal
+ :param Ry0: Magnetoresistance offset transverse
+ :param AMR_X: Anisotropic magnetoresistance longitudinal
+ :param AMR_Y: Anisotropic magnetoresistance transverse
+ :param SMR_X: Spin magnetoresistance longitudinal
+ :param SMR_Y: Spin magnetoresistance transverse
+ :param AHE: Anomalous Hall effect resistance offset (transverse only)
"""
@overload
- def __init__(*args, **kwargs) -> Any:
- ...
-
+ def __init__(*args, **kwargs) -> Any: ...
def clearLog(self) -> Dict[str, Any]:
"""
Reset current simulation state`
"""
...
- def getLayerMagnetisation(self, layer_id: str) -> CVector:
- ...
-
+ def getLayerMagnetisation(self, layer_id: str) -> CVector: ...
def getLog(self) -> Dict[str, List[float]]:
"""
Retrieve the simulation log [data].
"""
...
- def getMagnetoresistance(self) -> List[float]:
- ...
-
+ def getMagnetoresistance(self) -> List[float]: ...
def runSimulation(
self,
totalTime: float,
@@ -200,8 +253,9 @@ class Junction:
"""
...
- def setIECDriver(self, bottom_layer: str, top_layer: str,
- driver: ScalarDriver) -> None:
+ def setIECDriver(
+ self, bottom_layer: str, top_layer: str, driver: ScalarDriver
+ ) -> None:
"""Set IEC interaction between two layers.
The names of the params are only for convention. The IEC will be set
between bottomLyaer or topLayer, order is irrelevant.
@@ -210,8 +264,9 @@ class Junction:
"""
...
- def setQuadIECDriver(self, bottom_layer: str, top_layer: str,
- driver: ScalarDriver) -> None:
+ def setQuadIECDriver(
+ self, bottom_layer: str, top_layer: str, driver: ScalarDriver
+ ) -> None:
"""Set secondary (biquadratic term) IEC interaction between two layers.
The names of the params are only for convention. The IEC will be set
between bottomLyaer or topLayer, order is irrelevant.
@@ -220,48 +275,40 @@ class Junction:
"""
...
- def setLayerTemperatureDriver(self, layer_id: str,
- driver: ScalarDriver) -> None:
- ...
-
- def setLayerAnisotropyDriver(self, layer_id: str,
- driver: ScalarDriver) -> None:
- ...
-
- def setLayerCurrentDriver(self, layer_id: str,
- driver: ScalarDriver) -> None:
- ...
-
- def setLayerExternalFieldDriver(self, layer_id: str,
- driver: AxialDriver) -> None:
- ...
-
- def setLayerMagnetisation(self, layer_id: str, mag: CVector) -> None:
- ...
-
+ def setLayerTemperatureDriver(
+ self, layer_id: str, driver: ScalarDriver
+ ) -> None: ...
+ def setLayerAnisotropyDriver(self, layer_id: str, driver: ScalarDriver) -> None: ...
+ def setLayerCurrentDriver(self, layer_id: str, driver: ScalarDriver) -> None: ...
+ def setLayerExternalFieldDriver(
+ self, layer_id: str, driver: AxialDriver
+ ) -> None: ...
+ def setLayerMagnetisation(self, layer_id: str, mag: CVector) -> None: ...
@overload
- def setLayerOerstedFieldDriver(self, layer_id: str,
- driver: AxialDriver) -> None:
- ...
-
- def setLayerDampingLikeTorqueDriver(self, layer_id: str,
- driver: ScalarDriver) -> None:
+ def setLayerOerstedFieldDriver(
+ self, layer_id: str, driver: AxialDriver
+ ) -> None: ...
+ def setLayerDampingLikeTorqueDriver(
+ self, layer_id: str, driver: ScalarDriver
+ ) -> None:
"""Set the damping like torque driver for a layer.
:param layer_id: the layer id
:param driver: the driver
"""
...
- def setLayerFieldLikeTorqueDriver(self, layer_id: str,
- driver: ScalarDriver) -> None:
+ def setLayerFieldLikeTorqueDriver(
+ self, layer_id: str, driver: ScalarDriver
+ ) -> None:
"""Set the field like torque driver for a layer.
:param layer_id: the layer id
:param driver: the driver
"""
...
- def setLayerOneFNoise(self, layer_id: str, sources: int, bias: float,
- scale: float) -> None:
+ def setLayerOneFNoise(
+ self, layer_id: str, sources: int, bias: float, scale: float
+ ) -> None:
"""Set 1/f noise for a layer.
:param layer_id: the layer id
:param sources: the number of generation sources (the more the slower, but more acc.)
@@ -270,7 +317,6 @@ class Junction:
"""
...
-
class Layer:
def __init__(
self,
@@ -369,30 +415,20 @@ class Layer:
Automatically changes the solver to Euler-Heun."""
...
- def setExternalFieldDriver(self, driver: AxialDriver) -> None:
- ...
-
- def setMagnetisation(self, mag: CVector) -> None:
- ...
-
- def setOerstedFieldDriver(self, driver: AxialDriver) -> None:
- ...
-
- def setDampingLikeTorqueDriver(self, arg0: ScalarDriver) -> None:
+ def setExternalFieldDriver(self, driver: AxialDriver) -> None: ...
+ def setMagnetisation(self, mag: CVector) -> None: ...
+ def setOerstedFieldDriver(self, driver: AxialDriver) -> None: ...
+ def setDampingLikeTorqueDriver(self, driver: ScalarDriver) -> None:
"""Set a driver for the damping like torque of the layer."""
...
- def setFieldLikeTorqueDriver(self, arg0: ScalarDriver) -> None:
+ def setFieldLikeTorqueDriver(self, driver: ScalarDriver) -> None:
"""Set a driver for the field like torque of the layer."""
...
- def setReferenceLayer(self, ref: CVector) -> None:
- ...
-
+ def setReferenceLayer(self, ref: CVector) -> None: ...
@overload
- def setReferenceLayer(self, ref: "Reference") -> None:
- ...
-
+ def setReferenceLayer(self, ref: "Reference") -> None: ...
def setTopDipoleTensor(self, tensor: List[CVector]) -> None:
"""Set a dipole tensor from the top layer."""
...
@@ -405,12 +441,18 @@ class Layer:
"""Get Id of the layer"""
...
- def setAlternativeSTT(self, arg0: bool) -> None:
+ def setAlternativeSTT(self, setAlternative: bool) -> None:
"""Switch to an alternative STT forumulation (Taniguchi et al.)
https://iopscience.iop.org/article/10.7567/APEX.11.013005
"""
...
+ def setKappa(self, kappa: float) -> None:
+ """Set the kappa parameter for the layer -- determines SOT mixing
+ Hdl * kappa + Hfl
+ Allows you to turn off Hdl. Turning Hfl is via beta parameter.
+ """
+ ...
class NullDriver(ScalarDriver):
def __init__(self) -> None:
@@ -420,11 +462,8 @@ class NullDriver(ScalarDriver):
"""
...
-
class ScalarDriver:
- def __init__(self, *args, **kwargs) -> None:
- ...
-
+ def __init__(self, *args, **kwargs) -> None: ...
@staticmethod
def getConstantDriver(constantValue: float) -> "ScalarDriver":
"""
@@ -435,8 +474,9 @@ class ScalarDriver:
...
@staticmethod
- def getPulseDriver(constantValue: float, amplitude: "ScalarDriver",
- period: float, cycle: float) -> Any:
+ def getPulseDriver(
+ constantValue: float, amplitude: "ScalarDriver", period: float, cycle: float
+ ) -> Any:
"""
Produces a square pulse of certain period and cycle
:param constantValue: offset (vertical) of the pulse. The pulse amplitude will be added to this.
@@ -448,8 +488,9 @@ class ScalarDriver:
...
@staticmethod
- def getSineDriver(constantValue: float, amplitude: "ScalarDriver",
- frequency: float, phase: float) -> Any:
+ def getSineDriver(
+ constantValue: float, amplitude: "ScalarDriver", frequency: float, phase: float
+ ) -> Any:
"""
Produces a sinusoidal signal with some offset (constantValue), amplitude frequency and phase offset.
:param constantValue: vertical offset. The sine will oscillate around this value.
@@ -461,8 +502,9 @@ class ScalarDriver:
...
@staticmethod
- def getStepDriver(constantValue: float, amplitude: float, timeStart: float,
- timeStop: float) -> "ScalarDriver":
+ def getStepDriver(
+ constantValue: float, amplitude: float, timeStart: float, timeStop: float
+ ) -> ScalarDriver:
"""
Get a step driver. It has amplitude between timeStart and timeStop and 0 elsewhere
:param constantValue: offset of the pulse (vertical)
@@ -472,9 +514,14 @@ class ScalarDriver:
"""
...
- def getTrapezoidDriver(self, constantValue: float, amplitude: float,
- timeStart, edgeTime: float,
- steadyTime: float) -> Any:
+ @staticmethod
+ def getTrapezoidDriver(
+ constantValue: float,
+ amplitude: float,
+ timeStart,
+ edgeTime: float,
+ steadyTime: float,
+ ) -> ScalarDriver:
"""Create Trapezoid driver. Has a rising and a falling edge.
:param constantValue: offset of the pulse (vertical)
:param amplitude: amplitude that is added on top of the constantValue
@@ -484,22 +531,79 @@ class ScalarDriver:
"""
...
+ @staticmethod
+ def getGaussianImpulseDriver(
+ constantValue: float, amplitude: float, t0: float, sigma: float
+ ) -> ScalarDriver:
+ """Gaussian impulse driver. It has amplitude starts at t0 and falls off with sigma.
+
+ Formula:
+ A * exp(-((t - t0) ** 2) / (2 * sigma ** 2))
+
+ :param constantValue: offset of the pulse (vertical)
+ :param amplitude: amplitude that is added on top of the constantValue
+ :param t0: start of the pulse
+ :param sigma: fall-off of the Gaussian pulse
+ """
+ ...
+
+ @staticmethod
+ def getGaussianStepDriver(
+ constantValue: float, amplitude: float, t0: float, sigma: float
+ ) -> ScalarDriver:
+ """Gaussian step driver (erf function). It has amplitude starts at t0 and falls off with sigma.
+
+ Formula:
+ f(t) = constantValue + amplitude * (1 + erf((t - t0) / (sigma * sqrt(2))))
+
+ :param constantValue: offset of the pulse (vertical)
+ :param amplitude: amplitude that is added on top of the constantValue
+ :param t0: start of the pulse
+ :param sigma: fall-off of the Gaussian pulse
+ """
+ ...
+
+ @staticmethod
+ def getPosSineDriver(
+ constantValue: float, amplitude: float, frequency: float, phase: float
+ ) -> ScalarDriver:
+ """Produces a positive sinusoidal signal with some offset (constantValue), amplitude frequency and phase offset.
+ :param constantValue: vertical offset. The sine will oscillate around this value.
+ :param amplitude: amplitude of the sine wave
+ :param frequency: frequency of the sine
+ :param phase: phase of the sine in radians.
+ """
+ ...
+
+ @staticmethod
+ def getPulseDriver(
+ constantValue: float, amplitude: float, period: float, cycle: float
+ ) -> ScalarDriver:
+ """
+ Produces a square pulse of certain period and cycle
+ :param constantValue: offset (vertical) of the pulse. The pulse amplitude will be added to this.
+ :param amplitude: amplitude of the pulse signal
+ :param period: period of the signal in seconds
+ :param cycle: duty cycle of the signal -- a fraction between [0 and 1].
+ """
+ ...
class SolverMode:
"""SolverMode Indicator"""
+
DormandPrice: ClassVar[SolverMode] = ...
EulerHeun: ClassVar[SolverMode] = ...
RK4: ClassVar[SolverMode] = ...
-
+ Heun: ClassVar[SolverMode] = ...
class Reference:
"""Reference layer indicator."""
+
bottom: ClassVar[Reference] = ...
fixed: ClassVar[Reference] = ...
none: ClassVar[Reference] = ...
top: ClassVar[Reference] = ...
-
DormandPrice: SolverMode
EulerHeun: SolverMode
Heun: SolverMode
diff --git a/cmtj/models/domain_dynamics.py b/cmtj/models/domain_dynamics.py
index 1520f7e..5870c3e 100644
--- a/cmtj/models/domain_dynamics.py
+++ b/cmtj/models/domain_dynamics.py
@@ -1,5 +1,4 @@
import math
-from abc import ABC
from collections import defaultdict
from dataclasses import dataclass, field
from typing import Callable, List, Literal
diff --git a/cmtj/models/general_sb.py b/cmtj/models/general_sb.py
index 2ab4772..e4d4998 100644
--- a/cmtj/models/general_sb.py
+++ b/cmtj/models/general_sb.py
@@ -13,7 +13,7 @@
from ..utils import VectorObj, gamma, gamma_rad, mu0, perturb_position
from ..utils.solvers import RootFinder
-EPS = np.finfo('float64').resolution
+EPS = np.finfo("float64").resolution
def real_deocrator(fn):
@@ -47,8 +47,10 @@ def general_hessian_functional(N: int):
all_symbols.extend(
(sym.Symbol(r"\theta_" + indx_i), sym.Symbol(r"\phi_" + indx_i)))
energy_functional_expr = sym.Function("E")(*all_symbols)
- return get_hessian_from_energy_expr(
- N, energy_functional_expr), energy_functional_expr
+ return (
+ get_hessian_from_energy_expr(N, energy_functional_expr),
+ energy_functional_expr,
+ )
@lru_cache
@@ -67,10 +69,9 @@ def get_hessian_from_energy_expr(N: int, energy_functional_expr: sym.Expr):
indx_i = str(i)
# z = sym.Symbol("Z")
# these here must match the Ms symbols!
- z = sym.Symbol(
- r"\omega") * sym.Symbol(r"M_{" + indx_i + "}") * sym.sin(
- sym.Symbol(r"\theta_" + indx_i)) * sym.Symbol(r"t_{" + indx_i +
- "}")
+ z = (sym.Symbol(r"\omega") * sym.Symbol(r"M_{" + indx_i + "}") *
+ sym.sin(sym.Symbol(r"\theta_" + indx_i)) *
+ sym.Symbol(r"t_{" + indx_i + "}"))
for j in range(i, N):
# indx_j = str(j + 1) # for display purposes
indx_j = str(j)
@@ -127,11 +128,13 @@ def find_analytical_roots(N: int):
def get_all_second_derivatives(energy_functional_expr,
energy_expression,
- subs={}):
+ subs=None):
"""Get all second derivatives of the energy expression.
:param energy_functional_expr: symbolic energy_functional expression
:param energy_expression: symbolic energy expression (from solver)
:param subs: substitutions to be made."""
+ if subs is None:
+ subs = {}
second_derivatives = subs
symbols = energy_expression.free_symbols
for i, s1 in enumerate(symbols):
@@ -155,6 +158,7 @@ class LayerSB:
:param Ks: surface anisotropy (out-of plane, or perpendicular) value [J/m^3].
:param Ms: magnetisation saturation value in [A/m].
"""
+
_id: int
thickness: float
Kv: VectorObj
@@ -169,7 +173,7 @@ def __post_init__(self):
self.m = sym.ImmutableMatrix([
sym.sin(self.theta) * sym.cos(self.phi),
sym.sin(self.theta) * sym.sin(self.phi),
- sym.cos(self.theta)
+ sym.cos(self.theta),
])
def get_coord_sym(self):
@@ -181,9 +185,16 @@ def get_m_sym(self):
return self.m
@lru_cache(3)
- def symbolic_layer_energy(self, H: sym.ImmutableMatrix, J1top: float,
- J1bottom: float, J2top: float, J2bottom: float,
- top_layer: "LayerSB", down_layer: "LayerSB"):
+ def symbolic_layer_energy(
+ self,
+ H: sym.ImmutableMatrix,
+ J1top: float,
+ J1bottom: float,
+ J2top: float,
+ J2bottom: float,
+ top_layer: "LayerSB",
+ down_layer: "LayerSB",
+ ):
"""Returns the symbolic expression for the energy of the layer.
Coupling contribution comes only from the bottom layer (top-down crawl)"""
m = self.get_m_sym()
@@ -195,12 +206,13 @@ def symbolic_layer_energy(self, H: sym.ImmutableMatrix, J1top: float,
if top_layer is not None:
other_m = top_layer.get_m_sym()
- top_iec_energy = -(J1top / self.thickness) * m.dot(other_m) - (
- J2top / self.thickness) * m.dot(other_m)**2
+ top_iec_energy = (-(J1top / self.thickness) * m.dot(other_m) -
+ (J2top / self.thickness) * m.dot(other_m)**2)
if down_layer is not None:
other_m = down_layer.get_m_sym()
- bottom_iec_energy = -(J1bottom / self.thickness) * m.dot(
- other_m) - (J2bottom / self.thickness) * m.dot(other_m)**2
+ bottom_iec_energy = (
+ -(J1bottom / self.thickness) * m.dot(other_m) -
+ (J2bottom / self.thickness) * m.dot(other_m)**2)
return eng_non_interaction + top_iec_energy + bottom_iec_energy
def no_iec_symbolic_layer_energy(self, H: sym.ImmutableMatrix):
@@ -214,45 +226,56 @@ def no_iec_symbolic_layer_energy(self, H: sym.ImmutableMatrix):
field_energy = -mu0 * self.Ms * m.dot(H)
surface_anistropy = (-self.Ks +
- (1. / 2.) * mu0 * self.Ms**2) * (m[-1]**2)
+ (1.0 / 2.0) * mu0 * self.Ms**2) * (m[-1]**2)
volume_anisotropy = -self.Kv.mag * (m.dot(alpha)**2)
- return (field_energy + surface_anistropy + volume_anisotropy)
+ return field_energy + surface_anistropy + volume_anisotropy
def sb_correction(self):
- omega = sym.Symbol(r'\omega')
+ omega = sym.Symbol(r"\omega")
return (omega / gamma) * self.Ms * sym.sin(self.theta) * self.thickness
def __hash__(self) -> int:
return hash(str(self))
def __eq__(self, __value: "LayerSB") -> bool:
- return self._id == __value._id and self.thickness == __value.thickness and self.Kv == __value.Kv and self.Ks == __value.Ks and self.Ms == __value.Ms
+ return (self._id == __value._id and self.thickness == __value.thickness
+ and self.Kv == __value.Kv and self.Ks == __value.Ks
+ and self.Ms == __value.Ms)
@dataclass
class LayerDynamic(LayerSB):
alpha: float
- def rhs_llg(self, H: sym.Matrix, J1top: float, J1bottom: float,
- J2top: float, J2bottom: float, top_layer: "LayerSB",
- down_layer: "LayerSB"):
+ def rhs_llg(
+ self,
+ H: sym.Matrix,
+ J1top: float,
+ J1bottom: float,
+ J2top: float,
+ J2bottom: float,
+ top_layer: "LayerSB",
+ down_layer: "LayerSB",
+ ):
"""Returns the symbolic expression for the RHS of the spherical LLG equation.
Coupling contribution comes only from the bottom layer (top-down crawl)"""
- U = self.symbolic_layer_energy(H,
- J1top=J1top,
- J1bottom=J1bottom,
- J2top=J2top,
- J2bottom=J2bottom,
- top_layer=top_layer,
- down_layer=down_layer)
+ U = self.symbolic_layer_energy(
+ H,
+ J1top=J1top,
+ J1bottom=J1bottom,
+ J2top=J2top,
+ J2bottom=J2bottom,
+ top_layer=top_layer,
+ down_layer=down_layer,
+ )
# sum all components
- prefac = gamma_rad / (1. + self.alpha)**2
- inv_sin = 1. / (sym.sin(self.theta) + EPS)
+ prefac = gamma_rad / (1.0 + self.alpha)**2
+ inv_sin = 1.0 / (sym.sin(self.theta) + EPS)
dUdtheta = sym.diff(U, self.theta)
dUdphi = sym.diff(U, self.phi)
- dtheta = (-inv_sin * dUdphi - self.alpha * dUdtheta)
- dphi = (inv_sin * dUdtheta - self.alpha * dUdphi * (inv_sin)**2)
+ dtheta = -inv_sin * dUdphi - self.alpha * dUdtheta
+ dphi = inv_sin * dUdtheta - self.alpha * dUdphi * (inv_sin)**2
return prefac * sym.ImmutableMatrix([dtheta, dphi]) / self.Ms
def __eq__(self, __value: "LayerDynamic") -> bool:
@@ -264,10 +287,22 @@ def __hash__(self) -> int:
@dataclass
class Solver:
+ """General solver for the system.
+
+ :param layers: list of layers in the system.
+ :param J1: list of interlayer exchange constants. Goes (i)-(i+1), i = 0, 1, 2, ...
+ with i being the index of the layer.
+ :param J2: list of interlayer exchange constants.
+ :param H: external field.
+ :param Ndipole: list of dipole fields for each layer. Defaults to None.
+ Goes (i)-(i+1), i = 0, 1, 2, ... with i being the index of the layer.
+ """
+
layers: List[Union[LayerSB, LayerDynamic]]
J1: List[float]
J2: List[float]
H: VectorObj = None
+ Ndipole: List[List[VectorObj]] = None
def __post_init__(self):
if len(self.layers) != len(self.J1) + 1:
@@ -275,6 +310,19 @@ def __post_init__(self):
if len(self.layers) != len(self.J2) + 1:
raise ValueError("Number of layers must be 1 more than J2.")
+ self.dipoleMatrix: list[sym.Matrix] = None
+ if self.Ndipole is not None:
+ if len(self.layers) != len(self.Ndipole) + 1:
+ raise ValueError(
+ "Number of layers must be 1 more than number of tensors.")
+ if isinstance(self.layers[0], LayerDynamic):
+ raise ValueError(
+ "Dipole coupling is not yet supported for LayerDynamic.")
+ self.dipoleMatrix = [
+ sym.Matrix([d.get_cartesian() for d in dipole])
+ for dipole in self.Ndipole
+ ]
+
id_sets = {layer._id for layer in self.layers}
ideal_set = set(range(len(self.layers)))
if id_sets != ideal_set:
@@ -292,10 +340,12 @@ def get_layer_references(self, layer_indx, interaction_constant):
elif layer_indx == len(self.layers) - 1:
return self.layers[layer_indx -
1], None, interaction_constant[-1], 0
- return self.layers[layer_indx - 1], self.layers[
- layer_indx +
- 1], interaction_constant[layer_indx -
- 1], interaction_constant[layer_indx]
+ return (
+ self.layers[layer_indx - 1],
+ self.layers[layer_indx + 1],
+ interaction_constant[layer_indx - 1],
+ interaction_constant[layer_indx],
+ )
def compose_llg_jacobian(self, H: VectorObj):
"""Create a symbolic jacobian of the LLG equation in spherical coordinates."""
@@ -344,10 +394,15 @@ def create_energy(self,
if bottom_layer:
ratio_bottom = bottom_layer.thickness / (
layer.thickness + bottom_layer.thickness)
- energy += layer.symbolic_layer_energy(H, Jtop * ratio_top,
- Jbottom * ratio_bottom,
- J2top, J2bottom,
- top_layer, bottom_layer)
+ energy += layer.symbolic_layer_energy(
+ H,
+ Jtop * ratio_top,
+ Jbottom * ratio_bottom,
+ J2top,
+ J2bottom,
+ top_layer,
+ bottom_layer,
+ )
else:
# surface energy for correct angular gradient
for layer in self.layers:
@@ -358,10 +413,20 @@ def create_energy(self,
for i in range(len(self.layers) - 1):
l1m = self.layers[i].get_m_sym()
l2m = self.layers[i + 1].get_m_sym()
- ldot = (l1m.dot(l2m))
+ ldot = l1m.dot(l2m)
energy -= self.J1[i] * ldot
energy -= self.J2[i] * (ldot)**2
+ # dipole fields
+ if self.dipoleMatrix is not None:
+ mat = self.dipoleMatrix[i]
+ # is positive, just like demag
+ energy += ((mu0 / 2.0) * l1m.dot(mat * l2m) *
+ self.layers[i].Ms * self.layers[i + 1].Ms *
+ self.layers[i].thickness)
+ energy += ((mu0 / 2.0) * l2m.dot(mat * l1m) *
+ self.layers[i].Ms * self.layers[i + 1].Ms *
+ self.layers[i + 1].thickness)
return energy
def create_energy_hessian(self, equilibrium_position: List[float]):
@@ -415,14 +480,16 @@ def get_gradient_expr(self, accel="math"):
symbols.extend((theta, phi))
return sym.lambdify(symbols, grad_vector, accel)
- def adam_gradient_descent(self,
- init_position: np.ndarray,
- max_steps: int,
- tol: float = 1e-8,
- learning_rate: float = 1e-4,
- first_momentum_decay: float = 0.9,
- second_momentum_decay: float = 0.999,
- perturbation: float = 1e-6):
+ def adam_gradient_descent(
+ self,
+ init_position: np.ndarray,
+ max_steps: int,
+ tol: float = 1e-8,
+ learning_rate: float = 1e-4,
+ first_momentum_decay: float = 0.9,
+ second_momentum_decay: float = 0.999,
+ perturbation: float = 1e-6,
+ ):
"""
A naive implementation of Adam gradient descent.
See: ADAM: A METHOD FOR STOCHASTIC OPTIMIZATION, Kingma et Ba, 2015
@@ -440,24 +507,24 @@ def adam_gradient_descent(self,
m = np.zeros_like(current_position) # first momentum
v = np.zeros_like(current_position) # second momentum
eps = 1e-12
-
+ # history = []
while True:
step += 1
grad = np.asarray(gradfn(*current_position))
- m = first_momentum_decay * m + (1. - first_momentum_decay) * grad
- v = second_momentum_decay * v + (1. -
+ m = first_momentum_decay * m + (1.0 - first_momentum_decay) * grad
+ v = second_momentum_decay * v + (1.0 -
second_momentum_decay) * grad**2
- m_hat = m / (1. - first_momentum_decay**step)
- v_hat = v / (1. - second_momentum_decay**step)
+ m_hat = m / (1.0 - first_momentum_decay**step)
+ v_hat = v / (1.0 - second_momentum_decay**step)
new_position = current_position - learning_rate * m_hat / (
np.sqrt(v_hat) + eps)
if step > max_steps:
break
- # if np.linalg.norm(current_position - new_position) < tol:
- # break
if fast_norm(current_position - new_position) < tol:
break
current_position = new_position
+ # history.append(current_position)
+ # return np.asarray(current_position), np.asarray(history)
return np.asarray(current_position)
def single_layer_resonance(self, layer_indx: int, eq_position: np.ndarray):
@@ -480,19 +547,22 @@ def single_layer_resonance(self, layer_indx: int, eq_position: np.ndarray):
fmr = np.sqrt(float(fmr)) * gamma_rad / (2 * np.pi)
return fmr
- def solve(self,
- init_position: np.ndarray,
- max_steps: int = 1e9,
- learning_rate: float = 1e-4,
- adam_tol: float = 1e-8,
- first_momentum_decay: float = 0.9,
- second_momentum_decay: float = 0.999,
- perturbation: float = 1e-3,
- ftol: float = 0.01e9,
- max_freq: float = 80e9,
- force_single_layer: bool = False,
- force_sb: bool = False):
+ def solve(
+ self,
+ init_position: np.ndarray,
+ max_steps: int = 1e9,
+ learning_rate: float = 1e-4,
+ adam_tol: float = 1e-8,
+ first_momentum_decay: float = 0.9,
+ second_momentum_decay: float = 0.999,
+ perturbation: float = 1e-3,
+ ftol: float = 0.01e9,
+ max_freq: float = 80e9,
+ force_single_layer: bool = False,
+ force_sb: bool = False,
+ ):
"""Solves the system.
+ For dynamic LayerDynamic, the return is different, check :return.
1. Computes the energy functional.
2. Computes the gradient of the energy functional.
3. Performs a gradient descent to find the equilibrium position.
@@ -525,7 +595,8 @@ def solve(self,
learning_rate=learning_rate,
first_momentum_decay=first_momentum_decay,
second_momentum_decay=second_momentum_decay,
- perturbation=perturbation)
+ perturbation=perturbation,
+ )
if not force_sb and isinstance(self.layers[0], LayerDynamic):
eigenvalues, eigenvectors = self.dynamic_layer_solve(eq)
return eq, eigenvalues / 1e9, eigenvectors
@@ -561,9 +632,9 @@ def num_solve(self,
hes = self.create_energy_hessian(eq)
omega = sym.Symbol(r"\omega")
if len(self.layers) <= 3:
- y = real_deocrator(njit(sym.lambdify(omega, hes, 'math')))
+ y = real_deocrator(njit(sym.lambdify(omega, hes, "math")))
else:
- y = real_deocrator(sym.lambdify(omega, hes, 'math'))
+ y = real_deocrator(sym.lambdify(omega, hes, "math"))
r = RootFinder(0, max_freq, step=ftol, xtol=1e-8, root_dtype="float16")
roots = r.find(y)
# convert to GHz
@@ -625,7 +696,7 @@ def analytical_field_scan(
learning_rate: float = 1e-4,
first_momentum_decay: float = 0.9,
second_momentum_decay: float = 0.999,
- disable_tqdm: bool = False
+ disable_tqdm: bool = False,
) -> Iterable[Tuple[List[float], List[float], VectorObj]]:
"""Performs a field scan using the analytical solutions.
:param Hrange: the range of fields to scan.
diff --git a/cmtj/models/noise.py b/cmtj/models/noise.py
new file mode 100644
index 0000000..fa1b073
--- /dev/null
+++ b/cmtj/models/noise.py
@@ -0,0 +1,307 @@
+from typing import Literal
+
+import matplotlib.animation as animation
+import matplotlib.pyplot as plt
+import numpy as np
+from tqdm import tqdm
+
+
+def noise_model(
+ N: int,
+ steps: int = 3e5,
+ thermal_noise_std: float = 1e-3,
+ background_thermal_noise_std: float = 1e-3,
+ amplitude: float = 1e-3,
+ enable_oscillations: bool = False,
+ volume_distribution: Literal["pareto", "uniform"] = "pareto",
+ volume_distribution_params: dict = None,
+ freq_distribution: Literal["uniform", "functional"] = "uniform",
+ freq_distribution_params: dict = None,
+ frequency_scale: float = 1000,
+ time_scale: float = 1e-9,
+ phase_std: float = np.pi / 12,
+ dims: int = 1,
+ seed: int = 42,
+ offset: int = 1,
+ save_vectors: bool = False,
+ verbose: bool = False,
+ N_background_scale: int = 10,
+):
+ """
+ Generate a basic noise model.
+
+ :param N: Number of elements in the noise model.
+ :type N: int
+ :param steps: Number of simulation steps, defaults to 3e5.
+ :type steps: int, optional
+ :param thermal_noise_std: Standard deviation of the thermal noise, defaults to 1e-3.
+ :type thermal_noise_std: float, optional
+ :param amplitude: Amplitude of the oscillating frequencies, defaults to 1e-3.
+ :type amplitude: float, optional
+ :param enable_oscillations: Enable oscillations, defaults to False.
+ :type enable_oscillations: bool, optional
+ :param volume_distribution: Volume distribution type, defaults to "pareto".
+ :type volume_distribution: str, optional
+ :param volume_distribution_params: Parameters for the volume distribution, defaults to None.
+ :type volume_distribution_params: dict, optional
+ :param freq_distribution: Frequency distribution type, defaults to "uniform".
+ :type freq_distribution: str, optional
+ :param freq_distribution_params: Parameters for the frequency distribution, defaults to None.
+ :type freq_distribution_params: dict, optional
+ :param frequency_scale: Frequency scale, defaults to 1000.
+ :type frequency_scale: float, optional
+ :param time_scale: Time scale, defaults to 1e-9.
+ :type time_scale: float, optional
+ :param phase_std: Standard deviation of the phase, defaults to np.pi / 12.
+ :type phase_std: float, optional
+ :param dims: Number of dimensions, defaults to 1.
+ :type dims: int, optional
+ :param seed: Random seed, defaults to 42.
+ :type seed: int, optional
+ :param offset: Offset for the simulation, defaults to 1.
+ :type offset: int, optional
+ :param save_vectors: Save the stepwise m vectors per each volume, defaults to False.
+ :type save_vectors: bool, optional
+ :param verbose: Enable verbose mode, defaults to False.
+ :type verbose: bool, optional
+ :param N_background_scale: Number of background domains for oscillations and
+ background noise, defaults to 10xN.
+ :type N_background_scale: int, optional
+ :return: A tuple containing the following elements:
+ - m_values (ndarray): Array of shape (steps, dims) representing the simulated values.
+ - volumes (ndarray): Array of shape (N, 1) representing the volumes.
+ - freqs (ndarray): Array of shape (N,) representing the frequencies.
+ - time_scale (float): The time scale used in the simulation.
+ :rtype: tuple
+ """
+ if volume_distribution_params is None:
+ volume_distribution_params = {"shape": 1, "scale": 0.05}
+ if freq_distribution_params is None:
+ freq_distribution_params = {"low": 1, "high": 5000}
+ rng = np.random.default_rng(seed=seed)
+ if volume_distribution == "pareto":
+ volumes = rng.gamma(**volume_distribution_params, size=N)
+ elif volume_distribution == "uniform":
+ volumes = rng.random(N)
+ volumes = volumes / np.sum(volumes)
+
+ steps = int(steps)
+ volumes = volumes[:, np.newaxis]
+ volumes = np.sort(volumes, axis=0)
+
+ if freq_distribution == "functional":
+ freqs = (frequency_scale / volumes.ravel()).astype(int)[::-1] + 1
+ elif freq_distribution == "uniform":
+ freqs = rng.integers(**freq_distribution_params, size=N)
+
+ Np = N_background_scale * N
+ freqs_osc = rng.uniform(10, 1e11, Np)
+ phases = np.zeros(Np)
+ if phase_std > 0:
+ phases = rng.random(Np) * phase_std
+
+ vector_values = rng.random((N, dims))
+ m_values = np.zeros((steps, dims))
+ f_counts = np.zeros_like(freqs)
+ vectors = []
+ triggers = 0
+
+ def _oscillations(i: int):
+ return amplitude * np.sin(2 * np.pi * freqs_osc * i * time_scale +
+ phases).reshape(-1, 1)
+
+ def _background_noise(i: int):
+
+ return rng.normal(0, background_thermal_noise_std, dims)
+
+ if enable_oscillations and background_thermal_noise_std > 0:
+ raise ValueError(
+ "Cannot have both oscillations and background thermal noise enabled."
+ " Either set enable oscillations = False or background thermal noise = 0."
+ )
+ if enable_oscillations:
+ osc_fn = _oscillations
+ elif background_thermal_noise_std > 0:
+ osc_fn = _background_noise
+ for i in tqdm(range(offset, steps + offset), total=steps):
+ osc_vals = osc_fn(i).sum()
+ freq_disturbance = rng.integers(-3, 3, N)
+ freq_mask = i % (freqs + freq_disturbance) == 0
+ fsum = np.sum(freq_mask)
+ f_counts[freq_mask] += 1
+ if fsum > 0:
+ triggers += 1
+ vector_values[freq_mask] = rng.normal(0, thermal_noise_std,
+ (fsum, dims))
+ m_values[i - offset] += np.sum(volumes * vector_values,
+ axis=0) + osc_vals
+ else:
+ m_values[i - offset] = osc_vals + m_values[i - offset - 1]
+
+ if save_vectors:
+ vectors.append(vector_values.copy())
+
+ if verbose:
+ print(f"Triggers {triggers} out of {steps} steps")
+ if save_vectors:
+ return m_values, volumes, freqs, time_scale, f_counts, vectors
+ return m_values, volumes, freqs, time_scale, f_counts
+
+
+def autocorrelation(x, dT):
+ """
+ Compute the autocorrelation of the signal, based on the properties of the
+ power spectral density of the signal.
+ Taken from the StackOverflow answer:
+ https://stackoverflow.com/questions/643699/how-can-i-use-numpy-correlate-to-do-autocorrelation
+ :param x: the signal
+ :param dT: the time step
+ """
+ xp = x - np.mean(x)
+ f = np.fft.fft(xp)
+ p = np.abs(f)**2
+ pi = np.fft.ifft(p)
+ autocorr = np.real(pi)[:x.size // 2] / np.sum(xp**2)
+
+ # Create a lag array
+ lag = np.arange(0, len(autocorr)) * dT
+
+ return lag, autocorr
+
+
+def plot_noise_data(m_values: np.ndarray, volumes: np.ndarray,
+ freqs: np.ndarray, time_scale: float):
+ """
+ Plot noise data:
+ - Autocorrelation
+ - Volume vs Frequency
+ - Frequency vs Power
+ - Time vs Amplitude
+
+ :param m_values: An array of magnetic values.
+ :type m_values: np.ndarray
+ :param volumes: An array of volumes.
+ :type volumes: np.ndarray
+ :param freqs: An array of frequencies.
+ :type freqs: np.ndarray
+ :param time_scale: The time scale.
+ :type time_scale: float
+
+ :returns: The generated figure.
+ :rtype: matplotlib.figure.Figure
+ """
+ dims = m_values.shape[1]
+ R = m_values.dot(np.array([1, 0, 0])) if dims == 3 else m_values
+ lag, autocorr = autocorrelation(R.ravel(), time_scale)
+ k = m_values.shape[0]
+ with plt.style.context(["science", "nature"]):
+ w, h = plt.figaspect(1.0 / 4.0)
+ fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, dpi=300, figsize=(w, h))
+ ax1.plot(lag, autocorr, color="royalblue")
+ ax1.set_xlabel("Time lag (s)")
+ ax1.set_ylabel("Autocorrelation")
+ ax2.plot(volumes, freqs / 1000, color="crimson")
+ # histogram of volumes
+ ax25 = ax2.twinx()
+ ax25.hist(volumes,
+ bins=min(100, len(volumes)),
+ color="navy",
+ alpha=0.5,
+ label="Count")
+ ax25.set_ylabel("Count", rotation=-90, labelpad=10)
+ ax25.legend()
+ ax2.set_xlabel("Area (a.u.)")
+ ax2.set_ylabel("Modulo step activation (1000x)")
+ y = np.fft.fft(m_values, axis=0)
+ y = np.power(np.abs(y), 2)
+ y = y[:int(k // 2)]
+ x = np.fft.fftfreq(int(k), time_scale)
+ x = x[:int(k // 2)]
+ ax3.plot(x, y, color="royalblue")
+ ax3.set_xscale("log")
+ ax3.set_yscale("log")
+ ax3.set_xlabel("Frequency (Hz)")
+ ax3.set_ylabel("Power (a.u.)")
+ x_base = np.arange(0, k * time_scale, time_scale)
+ ax4.plot(x_base, m_values, color="forestgreen")
+ ax4.set_xlabel("Time")
+ ax4.set_ylabel("Amplitude")
+ fig.subplots_adjust(wspace=0.55)
+
+ # add letters
+ import matplotlib.transforms as mtransforms
+
+ for label, ax in zip("abcd", (ax1, ax2, ax3, ax4)):
+ # label physical distance in and down:
+ trans = mtransforms.ScaledTranslation(10 / 72, -5 / 72,
+ fig.dpi_scale_trans)
+ ax.text(
+ 0.0,
+ 1.0,
+ f"{label})",
+ transform=ax.transAxes + trans,
+ # fontsize="medium",
+ verticalalignment="top",
+ color="black",
+ bbox=dict(facecolor="none",
+ alpha=0.4,
+ edgecolor="none",
+ pad=3.0),
+ )
+
+ return fig
+
+
+def create_noise_animation(
+ volumes: list,
+ vector_values: np.ndarray,
+ save_type: str = "gif",
+ max_frames: int = 1000,
+):
+ """
+ Create a 2D noise animation given volumes/areas and vector values.
+
+ :param volumes: A list of volumes.
+ :type volumes: list
+ :param vector_values: An array of vector values.
+ :type vector_values: np.ndarray
+ :param save_type: The type of animation to save. Defaults to "gif".
+ :type save_type: str, optional
+ :param max_frames: The maximum number of frames. Defaults to 1000.
+ :type max_frames: int, optional
+ :returns: None
+ """
+
+ rng = np.random.default_rng(seed=42)
+ vector_values = np.asarray(vector_values).squeeze()
+ vector_values.shape
+ v = volumes.ravel()
+ v = v / v.sum()
+ n = 1000
+ xx, yy = np.meshgrid(np.arange(0, n), np.arange(0, n))
+ values = np.zeros_like(xx, dtype=np.float32)
+ volume_masks = []
+ for i, volume in enumerate(v):
+ x0, y0 = rng.integers(0, n, 2)
+ shape = (xx - x0)**2 + (yy - y0)**2
+ mask = shape <= (volume / np.pi) * ((n / 2)**2)
+ values[mask] = vector_values[105, i]
+ volume_masks.append(mask)
+
+ # create matplotlib animation
+ fig, ax = plt.subplots(dpi=300)
+ im = ax.imshow(values, cmap="viridis", vmin=-1, vmax=1)
+ # remove axes
+ ax.axis("off")
+
+ def update(step):
+ for i, mask in enumerate(volume_masks):
+ values[mask] = vector_values[step, i]
+ im.set_array(values)
+ ax.set_title(f"Step {step}")
+
+ ani = animation.FuncAnimation(fig, update, frames=max_frames, interval=0.8)
+ if save_type == "gif":
+ ani.save("animation.gif", writer="pillow")
+ else:
+ ani.save("animation.mp4", writer="ffmpeg")
diff --git a/cmtj/noise/__init__.pyi b/cmtj/noise/__init__.pyi
new file mode 100644
index 0000000..aa5c242
--- /dev/null
+++ b/cmtj/noise/__init__.pyi
@@ -0,0 +1,38 @@
+import cmtj
+
+class BufferedAlphaNoise:
+ """Create a buffer of alpha noise generator. Alpha can be in [0, 2]."""
+
+ def __init__(
+ self, bufferSize: int, alpha: float, std: float, scale: float
+ ) -> None: ...
+ def fillBuffer(self) -> None:
+ """Fill the buffer with the noise. This method is called only once."""
+ ...
+ def tick(self) -> float:
+ """Produce the next sample of the noise."""
+ ...
+ pass
+
+class VectorAlphaNoise:
+ """Create a vector alpha noise generator. Alpha can be in [0, 2]."""
+
+ def __init__(
+ self,
+ bufferSize: int,
+ alpha: float,
+ std: float,
+ scale: float,
+ axis: cmtj.Axis = cmtj.Axis.all,
+ ) -> None: ...
+ def getPrevSample(self) -> cmtj.CVector:
+ """Get the previous sample of the noise in a vector form."""
+ ...
+ def getScale(self) -> float:
+ """Get the scale of the noise."""
+ ...
+ def tick(self) -> float: ...
+ def tickVector(self) -> cmtj.CVector:
+ """Get the next sample of the noise in a vector form."""
+ ...
+ pass
diff --git a/cmtj/stack/__init__.pyi b/cmtj/stack/__init__.pyi
index b897528..5763093 100644
--- a/cmtj/stack/__init__.pyi
+++ b/cmtj/stack/__init__.pyi
@@ -71,6 +71,8 @@ class ParallelStack:
"""
...
+ def getMagnetisation(self, junction: int, layerId: str) -> cmtj.CVector: ...
+
class SeriesStack:
def __init__(self, junctionList: List[cmtj.Junction]) -> None:
"""
@@ -139,3 +141,4 @@ class SeriesStack:
:param mag: the magnetisation to be set.
"""
...
+ def getMagnetisation(self, junction: int, layerId: str) -> cmtj.CVector: ...
diff --git a/core/cvector.hpp b/core/cvector.hpp
index c576b6a..b7b5dd8 100644
--- a/core/cvector.hpp
+++ b/core/cvector.hpp
@@ -4,8 +4,19 @@
#include // for function
#include // for operator<<, ostream
#include // for runtime_error
-#include // for vector
-#include
+#include // for allocator, vector
+#include // for char_traits, basic_stringstream, basic_os..
+
+/// @brief A simple enum to represent the axis
+enum Axis
+{
+ xaxis,
+ yaxis,
+ zaxis,
+ all,
+ none
+};
+
template
class CVector
{
@@ -259,6 +270,13 @@ class CVector
ss << "[x:" << this->x << ", y:" << this->y << ", z:" << this->z << "]";
return ss.str();
}
+
+ const std::string toString() const
+ {
+ std::stringstream ss;
+ ss << "[x:" << this->x << ", y:" << this->y << ", z:" << this->z << "]";
+ return ss.str();
+ }
};
#endif // CORE_CVECTOR_HPP_
diff --git a/core/drivers.hpp b/core/drivers.hpp
index 657dbca..b9b3783 100644
--- a/core/drivers.hpp
+++ b/core/drivers.hpp
@@ -357,13 +357,6 @@ class NullDriver : public ScalarDriver
}
};
-enum Axis
-{
- xaxis,
- yaxis,
- zaxis
-};
-
template
class AxialDriver : public Driver
{
@@ -450,6 +443,10 @@ class AxialDriver : public Driver
return AxialDriver(NullDriver(), in, NullDriver());
case zaxis:
return AxialDriver(NullDriver(), NullDriver(), in);
+ case all:
+ return AxialDriver(in, in, in);
+ case none:
+ return AxialDriver(NullDriver(), NullDriver(), NullDriver());
}
return AxialDriver(NullDriver(), NullDriver(), NullDriver());
}
diff --git a/core/junction.hpp b/core/junction.hpp
index 4643c5b..44367b2 100644
--- a/core/junction.hpp
+++ b/core/junction.hpp
@@ -141,16 +141,11 @@ enum SolverMode
HEUN = 3
};
-// seems to be the faster so far.
-// static std::mt19937 generator((std::random_device{}()));
-
template
class Layer
{
private:
- OneFNoise* ofn;
-
ScalarDriver temperatureDriver;
ScalarDriver IECDriverTop;
ScalarDriver IECDriverBottom;
@@ -172,7 +167,7 @@ class Layer
Reference referenceType = NONE;
// the distribution is binded for faster generation
- // is also shared between 1f and Gaussian noise.
+ // is also shared between 1/f and Gaussian noise.
std::function distribution = std::bind(std::normal_distribution(0, 1), std::mt19937(std::random_device{}()));
CVector dWn, dWn2; // one for thermal, one for OneF
@@ -216,10 +211,21 @@ class Layer
dWn = CVector(this->distribution);
dWn.normalize();
this->cellVolume = this->cellSurface * this->thickness;
- this->ofn = new OneFNoise(0, 0., 0.);
+ this->ofn = std::shared_ptr>(new OneFNoise(0, 0., 0.));
}
public:
+ struct BufferedNoiseParameters
+ {
+ /* data */
+ T alphaNoise = 1.0;
+ T scaleNoise = 0.0;
+ T stdNoise = 0.0;
+ Axis axis = Axis::all;
+ };
+ BufferedNoiseParameters noiseParams;
+ std::shared_ptr> ofn;
+ std::shared_ptr> bfn;
bool includeSTT = false;
bool includeSOT = false;
@@ -440,9 +446,30 @@ class Layer
}
void setOneFNoise(unsigned int sources, T bias, T scale) {
- this->ofn = new OneFNoise(sources, bias, scale);
+ this->ofn = std::shared_ptr>(new OneFNoise(sources, bias, scale));
this->pinkNoiseSet = true;
- this->nonStochasticOneFSet = true; // by default turn it on, but in the stochastic sims, we will have to turn it off
+ // by default turn it on, but in the stochastic sims, we will have to turn it off
+ this->nonStochasticOneFSet = true;
+ }
+
+ void setAlphaNoise(T alpha, T std, T scale, Axis axis = Axis::all) {
+ if ((alpha < 0) || (alpha > 2))
+ throw std::runtime_error("alpha must be between 0 and 2");
+ this->noiseParams.alphaNoise = alpha;
+ this->noiseParams.stdNoise = std;
+ this->noiseParams.scaleNoise = scale;
+ this->noiseParams.axis = axis;
+ this->pinkNoiseSet = true;
+ }
+
+ void createBufferedAlphaNoise(unsigned int bufferSize) {
+ if (this->noiseParams.alphaNoise < 0)
+ throw std::runtime_error("alpha must be set before creating the noise!"
+ " Use setAlphaNoise function to set the alpha parameter.");
+
+ this->bfn = std::shared_ptr>(new VectorAlphaNoise(bufferSize,
+ this->noiseParams.alphaNoise,
+ this->noiseParams.stdNoise, this->noiseParams.scaleNoise, this->noiseParams.axis));
}
void setCurrentDriver(const ScalarDriver& driver)
@@ -678,6 +705,8 @@ class Layer
{
this->I_log = this->currentDriver.getCurrentScalarValue(time);
// use standard STT formulation
+ // see that literature reports Ms/MAGNETIC_PERMEABILITY
+ // but then the units don't match, we use Ms [T] which works
const T aJ = HBAR * this->I_log /
(ELECTRON_CHARGE * this->Ms * this->thickness);
// field like
@@ -724,6 +753,22 @@ class Layer
}
return dmdt * -GYRO * convTerm;
}
+
+ /**
+ * @brief Assumes the dW has the scale of sqrt(timeStep).
+ *
+ * @param currentMag
+ * @param dW - stochastic vector already scaled properly
+ * @return CVector
+ */
+ CVector stochasticTorque(const CVector& currentMag, const CVector& dW) {
+
+ const T convTerm = -GYRO / (1. + pow(this->damping, 2));
+ const CVector thcross = c_cross(currentMag, dW);
+ const CVector thcross2 = c_cross(currentMag, thcross);
+ return (thcross + thcross2 * this->damping) * convTerm;
+ }
+
const CVector calculateLLGWithFieldTorqueDipoleInjection(T time, const CVector& m,
const CVector& bottom, const CVector& top,
const CVector& dipole, T timeStep, const CVector& Hfluctuation = CVector())
@@ -799,65 +844,6 @@ class Layer
this->mag = m_t;
}
- void dormandPriceStep(T time, T timeStep, const CVector& bottom, const CVector& top)
- {
- CVector m_t = this->mag;
- CVector e_t;
- std::array, 7> K;
- // if (this->hopt < 0)
- // {
- // }
- // this makes the step non-adaptive
- // we will deal with this problem later, this requires consistency in step
- // across all the layers
- this->hopt = timeStep;
- // define Buchter tableau below
- // there are redundant zeros, but for clarity, we'll leave them
- const std::array c = {
- 0., 1. / 5., 3. / 10., 4. / 5., 8. / 9., 1., 1. };
- const std::array b = {
- 35. / 384., 0, 500. / 1113., 125. / 192., -2187. / 6784., 11. / 84., 0. };
- const std::array b2 = {
- 5179. / 57600., 0, 7571. / 16695., 393. / 640., -92097. / 339200., 187. / 2100., 1. / 40. };
- // extra braces are required even though the struct is 2D only
- const std::array, 7> aCoefs = { {{0, 0, 0, 0, 0, 0, 0},
- {1. / 5., 0, 0, 0, 0, 0, 0},
- {3. / 40., 9. / 40., 0, 0, 0, 0, 0},
- {44. / 45., -56. / 15., 32. / 9., 0, 0, 0, 0},
- {19372. / 6561., -25360. / 2187., 64448. / 6561., -212. / 729., 0, 0, 0},
- {9017. / 3168., -355. / 33., 46732. / 5247., 49. / 176., -5103. / 18656., 0, 0},
- {35. / 384., 0., 500. / 1113., 125. / 192., -2187. / 6784., 11. / 84., 0.}} };
- // compute the first
- K[0] = calculateLLGWithFieldTorque(time, this->mag, bottom, top, this->hopt);
- m_t = m_t + K[0] * b[0] * this->hopt;
- for (int i = 1; i < 7; i++)
- {
- CVector kS;
- for (int j = 0; j < i; j++)
- {
- kS = kS + K[j] * aCoefs[i][j];
- }
- K[i] = calculateLLGWithFieldTorque(time + c[i] * this->hopt, this->mag + kS * this->hopt,
- bottom, top, this->hopt);
- m_t = m_t + K[i] * b[i] * this->hopt; // this is function estimate
- e_t = e_t + K[i] * (b[i] - b2[i]) * this->hopt; // this is error estimate
- }
- // adapt the step size
- const T eps = 1e-3;
- const T s = pow((eps * timeStep) / (2 * e_t.length()), 1. / 5.);
- this->hopt = s * timeStep;
- m_t.normalize();
- this->mag = m_t;
- if (isnan(this->mag.x))
- {
- throw std::runtime_error("NAN magnetisation");
- }
- }
-
- CVector non_stochastic_llg(const CVector& cm, T time, T timeStep, const CVector& bottom, const CVector& top)
- {
- return calculateLLGWithFieldTorque(time, cm, bottom, top, timeStep);
- }
CVector stochastic_llg(const CVector& cm, T time, T timeStep,
const CVector& bottom, const CVector& top, const CVector& dW, const CVector& dW2, const T& HoneF)
@@ -880,6 +866,10 @@ class Layer
const T getStochasticOneFNoise(T time) {
if (!this->pinkNoiseSet)
return 0;
+ else if (this->noiseParams.scaleNoise != 0) {
+ // use buffered noise if available
+ return this->bfn->tick();
+ }
return this->ofn->tick();
}
@@ -888,112 +878,25 @@ class Layer
if (this->cellVolume == 0.0)
throw std::runtime_error("Cell surface cannot be 0 during temp. calculations!");
const T currentTemp = this->temperatureDriver.getCurrentScalarValue(time);
- const T mainFactor = (2 * this->damping * MAGNETIC_PERMEABILITY * BOLTZMANN_CONST * currentTemp) / (this->Ms * this->cellVolume * timeStep);
+ const T mainFactor = (2 * this->damping * BOLTZMANN_CONST * currentTemp) / (this->Ms * this->cellVolume * GYRO);
return sqrt(mainFactor);
}
- CVector nonStochasticLangevin(T time, T timeStep)
+ CVector getStochasticLangevinVector(const T& time, const T& timeStep)
{
+ if (!this->temperatureSet)
+ return CVector();
const T Hthermal_temp = this->getLangevinStochasticStandardDeviation(time, timeStep);
const CVector dW = CVector(this->distribution);
return dW * Hthermal_temp;
}
- CVector nonStochasticOneFNoise(T time, T timestep) {
- const T pinkNoise = this->ofn->tick();
- const CVector dW2 = CVector(this->distribution);
- return dW2 * pinkNoise;
- }
-
- /**
- * @brief Computes a single Euler-Heun step [DEPRECATED].
- * [DEPRECATED] This is the old Euler-Heun method, Heun is preferred.
- * Bottom and top are relative to the current layer.
- * They are used to compute interactions.
- * @param time: current time of the simulation
- * @param timeStep: integration time of the solver
- * @param bottom: bottom layer to the current layer
- * @param top: top layer to the current layer
- */
- void euler_heun_step(T time, T timeStep, const CVector& bottom, const CVector& top)
- {
- // we compute the two below in stochastic part, not non stochastic.
- this->nonStochasticTempSet = false;
- this->nonStochasticOneFSet = false;
- // this is Stratonovich integral
- if (isnan(this->mag.x))
- {
- throw std::runtime_error("NAN magnetisation");
- }
- // Brownian motion sample
- // Generate the noise from the Brownian motion
- // dW2 is used for 1/f noise generation
- CVector dW = CVector(this->distribution); // * sqrt(timeStep);
- CVector dW2 = CVector(this->distribution); //* sqrt(timeStep);
- // squared dW -- just utility
- dW.normalize();
- dW2.normalize();
- // f_n is the vector of non-stochastic part at step n
- // multiply by timeStep (h) here for convenience
- const T Honef = this->getStochasticOneFNoise(time);
- const CVector f_n = non_stochastic_llg(this->mag, time, timeStep, bottom, top) * timeStep;
- // g_n is the stochastic part of the LLG at step n
- const CVector g_n = stochastic_llg(this->mag, time, timeStep, bottom, top, dW, dW2, Honef) * timeStep;
-
- // actual solution
- // approximate next step ytilde
- const CVector mapprox = this->mag + g_n;
- // calculate the approx g_n
- const CVector g_n_approx = stochastic_llg(mapprox, time, timeStep, bottom, top, dW, dW2, Honef) * timeStep;
- // CVector m_t = this->mag + f_n + g_n + (g_n_approx - g_n) * 0.5;
- CVector m_t = this->mag + f_n + (g_n_approx + g_n) * 0.5;
- m_t.normalize();
- this->mag = m_t;
- }
-
- /**
- * @brief Computes a single Heun step.
- * This method is preferred over Euler-Heun method.
- * Bottom and top are relative to the current layer.
- * They are used to compute interactions.
- * @param time: current time of the simulation
- * @param timeStep: integration time of the solver
- * @param bottom: bottom layer to the current layer
- * @param top: top layer to the current layer
- */
- void heun_step(T time, T timeStep, const CVector& bottom, const CVector& top) {
- // we compute the two below in stochastic part, not non stochastic.
- this->nonStochasticTempSet = false;
- this->nonStochasticOneFSet = false;
- // this is Stratonovich integral
- if (isnan(this->mag.x))
- {
- throw std::runtime_error("NAN magnetisation");
+ CVector getOneFVector() {
+ if (this->noiseParams.scaleNoise != 0) {
+ // use buffered noise if available
+ return this->bfn->tickVector();
}
- // Brownian motion sample
- // Generate the noise from the Brownian motion
- const T Honef_scale = this->getStochasticOneFNoise(time);
- const T Hthermal_scale = this->getLangevinStochasticStandardDeviation(time, timeStep);
- const CVector Hlangevin = dWn * Hthermal_scale;
- const CVector Honef = dWn2 * Honef_scale;
- const CVector m_t = this->mag;
- const CVector f_n = this->calculateLLGWithFieldTorque(time, m_t, bottom, top, timeStep, Hlangevin + Honef);
- // immediate m approximation
- CVector m_approx = m_t + f_n * timeStep;
- CVector dW = CVector(this->distribution);
- CVector dW2 = CVector(this->distribution);
- dW.normalize();
- dW2.normalize();
- m_approx.normalize();
- const CVector Hlangevin_approx = dW * Hthermal_scale;
- const CVector Honef_approx = dW2 * Honef_scale;
- const CVector f_approx = this->calculateLLGWithFieldTorque(time + timeStep,
- m_approx, bottom, top, timeStep, Hlangevin_approx + Honef_approx);
- dWn = dW; // replace
- dWn2 = dW2;
- CVector nm_t = this->mag + (f_n + f_approx) * 0.5 * timeStep;
- nm_t.normalize();
- this->mag = nm_t;
+ return CVector();
}
};
@@ -1017,7 +920,7 @@ class Junction
std::vector Rx0, Ry0, AMR_X, AMR_Y, SMR_X, SMR_Y, AHE;
std::unordered_map> log;
- // std::string fileSave;
+
unsigned int logLength = 0;
unsigned int layerNo;
std::string Rtag = "R";
@@ -1060,7 +963,8 @@ class Junction
this->Rap = Rap;
this->MR_mode = CLASSIC;
// A string representing the tag for the junction's resistance value.
- this->Rtag = "R_" + this->layers[0].id + "_" + this->layers[1].id;
+ if (this->layerNo == 2)
+ this->Rtag = "R_" + this->layers[0].id + "_" + this->layers[1].id;
}
/**
@@ -1113,32 +1017,6 @@ class Junction
this->MR_mode = STRIP;
}
- /**
- * @brief Select a solver for the LLGS/sLLGS solutions
- * Available: DormandPrice, EulerHeun, RK4
- * @param solverMode one of EULER_HEUN, DORMAND_PRICE or RK4
- */
- const auto selectSolver(SolverMode solverMode)
- {
- switch (solverMode)
- {
- case EULER_HEUN:
- return &Layer::euler_heun_step;
-
- case DORMAND_PRICE:
- return &Layer::dormandPriceStep;
-
- case RK4:
- return &Layer::rk4_step;
-
- case HEUN:
- return &Layer::heun_step;
-
- default:
- return &Layer::rk4_step;
- }
- }
-
/**
* @brief Get Ids of the layers in the junction.
* @return vector of layer ids.
@@ -1180,7 +1058,7 @@ class Junction
}
if (!found)
{
- throw std::runtime_error("Failed to find a layer with a given id!");
+ throw std::runtime_error("Failed to find a layer with a given id: " + layerID + "!");
}
}
void axiallayerSetter(const std::string& layerID, axialDriverSetter functor, AxialDriver driver)
@@ -1196,7 +1074,7 @@ class Junction
}
if (!found)
{
- throw std::runtime_error("Failed to find a layer with a given id!");
+ throw std::runtime_error("Failed to find a layer with a given id: " + layerID + "!");
}
}
void setLayerTemperatureDriver(const std::string& layerID, const ScalarDriver& driver)
@@ -1288,7 +1166,7 @@ class Junction
}
if (!found)
{
- throw std::runtime_error("Failed to match the layer order or find layer ids!");
+ throw std::runtime_error("Failed to match the layer order or find layer ids: " + bottomLayer + " and " + topLayer + "!");
}
}
@@ -1315,7 +1193,7 @@ class Junction
}
if (!found)
{
- throw std::runtime_error("Failed to match the layer order or find layer ids!");
+ throw std::runtime_error("Failed to match the layer order or find layer ids: " + bottomLayer + " and " + topLayer + "!");
}
}
@@ -1332,7 +1210,7 @@ class Junction
}
if (!found)
{
- throw std::runtime_error("Failed to find a layer with a given id!");
+ throw std::runtime_error("Failed to find a layer with a given id: " + layerID + "!");
}
}
@@ -1381,7 +1259,7 @@ class Junction
if (res != this->layers.end()) {
return *res;
}
- throw std::runtime_error("Failed to find a layer with a given id!");
+ throw std::runtime_error("Failed to find a layer with a given id " + layerID + "!");
}
/**
@@ -1482,7 +1360,7 @@ class Junction
}
typedef void (Layer::* solverFn)(T t, T timeStep, const CVector& bottom, const CVector& top);
-
+ typedef void (Junction::* runnerFn)(solverFn& functor, T& t, T& timeStep);
/**
* @brief Run Euler-Heun or RK4 method for a single layer.
*
@@ -1525,6 +1403,101 @@ class Junction
}
}
+ void eulerHeunSolverStep(solverFn& functor, T& t, T& timeStep) {
+ /*
+ Euler Heun method (stochastic heun)
+
+ y_np = y + g(y,t,dW)*dt
+ g_sp = g(y_np,t+1,dW)
+ y(t+1) = y + dt*f(y,t) + .5*(g(y,t,dW)+g_sp)*sqrt(dt)
+
+ with f being the non-stochastic part and g the stochastic part
+ */
+ // draw the noise for each layer, dW
+ std::vector> mPrime(this->layerNo, CVector());
+ for (unsigned int i = 0; i < this->layerNo; i++) {
+ // todo: after you're done, double check the thermal magnitude and dt scaling there
+ const CVector dW = this->layers[i].getStochasticLangevinVector(t, timeStep) + this->layers[i].getOneFVector();
+ const CVector bottom = (i == 0) ? CVector() : this->layers[i - 1].mag;
+ const CVector top = (i == this->layerNo - 1) ? CVector() : this->layers[i + 1].mag;
+
+ const CVector fnApprox = this->layers[i].calculateLLGWithFieldTorque(
+ t, this->layers[i].mag, bottom, top, timeStep);
+ const CVector gnApprox = this->layers[i].stochasticTorque(this->layers[i].mag, dW);
+
+ // theoretically we have 2 options
+ // 1. calculate only the stochastic part with the second approximation
+ // 2. calculate the second approximation of m with the stochastic and non-stochastic
+ // part and then use if for torque est.
+ const CVector mNext = this->layers[i].mag + gnApprox * sqrt(timeStep);
+ const CVector gnPrimeApprox = this->layers[i].stochasticTorque(mNext, dW);
+ mPrime[i] = this->layers[i].mag + fnApprox * timeStep + 0.5 * (gnApprox + gnPrimeApprox) * sqrt(timeStep);
+ }
+
+ for (unsigned int i = 0; i < this->layerNo; i++) {
+ this->layers[i].mag = mPrime[i];
+ this->layers[i].mag.normalize();
+ }
+ }
+
+
+ void heunSolverStep(solverFn& functor, T& t, T& timeStep) {
+ /*
+ Heun method
+ y'(t+1) = y(t) + dy(y, t)
+ y(t+1) = y(t) + 0.5 * (dy(y, t) + dy(y'(t+1), t+1))
+ */
+ /*
+ Stochastic Heun method
+ y_np = y + g(y,t,dW)*dt
+ g_sp = g(y_np,t+1,dW)
+ y' = y_n + f_n * dt + g_n * dt
+ f' = f(y, )
+ y(t+1) = y + dt*f(y,t) + .5*(g(y,t,dW)+g_sp)*sqrt(dt)
+ */
+ std::vector> fn(this->layerNo, CVector());
+ std::vector> gn(this->layerNo, CVector());
+ std::vector> dW(this->layerNo, CVector());
+ std::vector> mNext(this->layerNo, CVector());
+ // first approximation
+
+ // make sure that
+ // 1. Thermal field is added if needed
+ // 2. One/f noise is added if needed
+ // 3. The timestep is correctly multiplied
+
+ for (unsigned int i = 0; i < this->layerNo; i++)
+ {
+ const CVector bottom = (i == 0) ? CVector() : this->layers[i - 1].mag;
+ const CVector top = (i == this->layerNo - 1) ? CVector() : this->layers[i + 1].mag;
+
+ fn[i] = this->layers[i].calculateLLGWithFieldTorque(
+ t, this->layers[i].mag, bottom, top, timeStep);
+
+ // draw the noise for each layer, dW
+ dW[i] = this->layers[i].getStochasticLangevinVector(t, timeStep) + this->layers[i].getOneFVector();
+ gn[i] = this->layers[i].stochasticTorque(this->layers[i].mag, dW[i]);
+
+ mNext[i] = this->layers[i].mag + fn[i] * timeStep + gn[i] * sqrt(timeStep);
+ }
+ // second approximation
+ for (unsigned int i = 0; i < this->layerNo; i++)
+ {
+ const CVector bottom = (i == 0) ? CVector() : mNext[i - 1];
+ const CVector top = (i == this->layerNo - 1) ? CVector() : mNext[i + 1];
+ // first approximation is already multiplied by timeStep
+ this->layers[i].mag = this->layers[i].mag + 0.5 * timeStep * (
+ fn[i] + this->layers[i].calculateLLGWithFieldTorque(
+ t + timeStep, mNext[i],
+ bottom,
+ top, timeStep)
+ ) + 0.5 * (gn[i] + this->layers[i].stochasticTorque(mNext[i], dW[i])) * sqrt(timeStep);
+ // normalise
+ this->layers[i].mag.normalize();
+ }
+
+ }
+
/**
* @brief Calculate strip magnetoresistance for multilayer.
*
@@ -1606,6 +1579,45 @@ class Junction
}
}
+ std::tuple getSolver(SolverMode mode, unsigned int totalIterations) {
+ SolverMode localMode = mode;
+ for (auto& l : this->layers)
+ {
+ if (l.hasTemperature())
+ {
+ // if at least one temp. driver is set
+ // then use heun for consistency
+ if (localMode != HEUN && localMode != EULER_HEUN) {
+ std::cout << "[WARNING] Solver automatically changed to Euler Heun for stochastic calculation." << std::endl;
+ localMode = EULER_HEUN;
+ }
+ }
+ if (l.noiseParams.scaleNoise != 0) {
+ // if at least one temp. driver is set
+ // then use heun for consistency
+ if (localMode != HEUN && localMode != EULER_HEUN) {
+ std::cout << "[WARNING] Solver automatically changed to Euler Heun for stochastic calculation." << std::endl;
+ localMode = EULER_HEUN;
+ }
+ // create a buffer
+ l.createBufferedAlphaNoise(totalIterations);
+ }
+ }
+ auto solver = &Layer::rk4_step;
+
+ // assign a runner function pointer from junction
+ auto runner = &Junction::runMultiLayerSolver;
+ if (this->layerNo == 1)
+ runner = &Junction::runSingleLayerSolver;
+ if (localMode == HEUN)
+ runner = &Junction::heunSolverStep;
+ else if (localMode == EULER_HEUN)
+ runner = &Junction::eulerHeunSolverStep;
+
+ return std::make_tuple(runner, solver, localMode);
+ }
+
+
/**
* Main run simulation function. Use it to run the simulation.
* @param totalTime: total time of a simulation, give it in seconds. Typical length is in ~couple ns.
@@ -1628,33 +1640,13 @@ class Junction
const unsigned int totalIterations = (int)(totalTime / timeStep);
const unsigned int writeEvery = (int)(writeFrequency / timeStep);
std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
- auto solver = this->selectSolver(mode);
// pick a solver based on drivers
- for (auto& l : this->layers)
- {
- if (l.hasTemperature())
- {
- if (mode != HEUN && mode != EULER_HEUN) {
- std::cout << "[WARNING] Solver automatically changed to Heun for stochastic calculation." << std::endl;
- }
- // if at least one temp. driver is set
- // then use heun for consistency
- solver = this->selectSolver(HEUN);
- break;
- }
- }
+ auto [runner, solver, _] = getSolver(mode, totalIterations);
for (unsigned int i = 0; i < totalIterations; i++)
{
T t = i * timeStep;
- if (this->layerNo == 1)
- {
- runSingleLayerSolver(solver, t, timeStep);
- }
- else
- {
- runMultiLayerSolver(solver, t, timeStep);
- }
+ (*this.*runner)(solver, t, timeStep);
if (!(i % writeEvery))
{
diff --git a/core/llgb.hpp b/core/llgb.hpp
new file mode 100644
index 0000000..57e58d0
--- /dev/null
+++ b/core/llgb.hpp
@@ -0,0 +1,571 @@
+#include "junction.hpp"
+#include