diff --git a/.gitignore b/.gitignore index c0c5f57d..d6260fcb 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,11 @@ *.pdf *.tmp *.so +*.o +*.dSYM +fidimag/**/*.cpp + + local/ *~ *.bak diff --git a/Makefile b/Makefile index 3e597546..33a36c6f 100644 --- a/Makefile +++ b/Makefile @@ -3,17 +3,50 @@ EXTENSIONS_DIR = ${PROJECT_DIR}/fidimag/extensions PYTHON = python3 PYTEST = ${PYTHON} -m pytest -##################### -# Cython Extensions # -##################### +FFTW_INC=local/include +FFTW_LIB=local/lib + +SUNDIALS_INC=local/include +SUNDIALS_LIB=local/lib + +CPPFLAGS = -fPIC -g -Iinclude/ \ + -I${FFTW_INC} -L${FFTW_LIB} \ + -I${SUNDIALS_INC} -L${SUNDIALS_LIB} \ + -Inative/include \ + -lm \ + -lfftw3_omp -lfftw3 \ + -lsundials_cvodes -lsundials_nvecserial -lsundials_nvecopenmp \ + -lblas -llapack \ + -fopenmp \ + -std=c++14 + +LDFLAGS = -shared +SOURCES = $(shell echo native/src/*.cpp) +OBJECTS = $(SOURCES:.cpp=.o) +LIBRARY = local/lib/libfidimag.so + +######## +# Build +######## + + +all: $(LIBRARY) build + + +$(LIBRARY) : $(OBJECTS) + $(CXX) $(CPPFLAGS) $(OBJECTS) -o $@ $(LDFLAGS) + + +build: $(LIBRARY) + CC=${CC} CXX=${CXX} CPPFLAGS="${CPPFLAGS}" ${PYTHON} setup.py build_ext --inplace -build: - ${PYTHON} setup.py build_ext --inplace clean: rm -rf ${EXTENSIONS_DIR}/* touch ${EXTENSIONS_DIR}/__init__.py rm -rf build + rm -rf $(OBJECTS) $(TARGET) *.dSYM + docker: docker build -t fidimag -f ./docker/travis/Dockerfile . diff --git a/docker/travis/Dockerfile b/docker/travis/Dockerfile index f326a156..86cb9520 100644 --- a/docker/travis/Dockerfile +++ b/docker/travis/Dockerfile @@ -1,4 +1,4 @@ -FROM ubuntu:16.04 +FROM ubuntu:18.04 # To build this image `docker build -t fidimag .` # Then you can drop into a live bash shell with `docker run -it fidimag`. diff --git a/fidimag/__init__.py b/fidimag/__init__.py index c8aa8bdc..7e22183d 100644 --- a/fidimag/__init__.py +++ b/fidimag/__init__.py @@ -48,6 +48,3 @@ Processor: {} """ print(message.format(FIDIMAG_DIR, e, platform.machine(), platform.platform(), platform.processor())) - - -__version__ = '0.2' diff --git a/fidimag/atomistic/lib/clib.pyx b/fidimag/atomistic/a_clib.pyx similarity index 98% rename from fidimag/atomistic/lib/clib.pyx rename to fidimag/atomistic/a_clib.pyx index b523b2c2..032ec26a 100644 --- a/fidimag/atomistic/lib/clib.pyx +++ b/fidimag/atomistic/a_clib.pyx @@ -1,8 +1,12 @@ +# distutils: language = c++ import numpy cimport numpy as np +cimport fidimag.common.vectormath as vmath np.import_array() -cdef extern from "fidimag_random.h": + + +cdef extern from "a_random.h": ctypedef struct mt19937_state: pass @@ -18,7 +22,7 @@ cdef extern from "time.h": ctypedef int time_t time_t time(time_t *timer) -cdef extern from "clib.h": +cdef extern from "a_clib.h": void run_step_mc(mt19937_state *state, double *spin, double *new_spin, int *ngbs, int *nngbs, int n_ngbs, double J, double J1, double D, double D1, @@ -95,8 +99,6 @@ cdef extern from "clib.h": # ------------------------------------------------------------------------- - void normalise(double *m, int *pins, int n) - # used for sllg void llg_rhs_dw_c(double *m, double *h, double *dm, double *T, double *alpha, double *mu_s_inv, int *pins, double *eta, int n, @@ -270,7 +272,7 @@ def compute_demag_full(double [:] spin, def normalise_spin(np.ndarray[double, ndim=1, mode="c"] spin, np.ndarray[int, ndim=1, mode="c"] pins, n): - normalise(&spin[0], &pins[0], n) + vmath.normalise(&spin[0], &pins[0], n) def compute_llg_rhs_dw(double [:] dm, double [:] spin, diff --git a/fidimag/atomistic/anisotropy.py b/fidimag/atomistic/anisotropy.py index 35048487..3f7c34e3 100644 --- a/fidimag/atomistic/anisotropy.py +++ b/fidimag/atomistic/anisotropy.py @@ -1,4 +1,4 @@ -import fidimag.extensions.clib as clib +import fidimag.extensions.a_clib as clib import numpy as np from .energy import Energy import fidimag.common.helper as helper diff --git a/fidimag/atomistic/atomistic_driver.py b/fidimag/atomistic/atomistic_driver.py index f09ac49d..bcd39dcb 100644 --- a/fidimag/atomistic/atomistic_driver.py +++ b/fidimag/atomistic/atomistic_driver.py @@ -3,7 +3,7 @@ from fidimag.common.driver_base import DriverBase -import fidimag.extensions.clib as clib +import fidimag.extensions.a_clib as clib import fidimag.extensions.cvode as cvode from fidimag.common.vtk import VTK import fidimag.common.constant as const diff --git a/fidimag/atomistic/demag_full.py b/fidimag/atomistic/demag_full.py index f2db488b..eaee35a3 100644 --- a/fidimag/atomistic/demag_full.py +++ b/fidimag/atomistic/demag_full.py @@ -1,4 +1,4 @@ -import fidimag.extensions.clib as clib +import fidimag.extensions.a_clib as clib import numpy as np from .energy import Energy diff --git a/fidimag/atomistic/dmi.py b/fidimag/atomistic/dmi.py index d15f2aaf..617d07f6 100644 --- a/fidimag/atomistic/dmi.py +++ b/fidimag/atomistic/dmi.py @@ -1,4 +1,4 @@ -import fidimag.extensions.clib as clib +import fidimag.extensions.a_clib as clib from .energy import Energy import numpy as np diff --git a/fidimag/atomistic/exchange.py b/fidimag/atomistic/exchange.py index 4d4ed3ec..1a46c343 100644 --- a/fidimag/atomistic/exchange.py +++ b/fidimag/atomistic/exchange.py @@ -1,4 +1,4 @@ -import fidimag.extensions.clib as clib +import fidimag.extensions.a_clib as clib import numpy as np from .energy import Energy diff --git a/fidimag/atomistic/hexagonal_mesh.py b/fidimag/atomistic/hexagonal_mesh.py index ee6ad6da..77420938 100644 --- a/fidimag/atomistic/hexagonal_mesh.py +++ b/fidimag/atomistic/hexagonal_mesh.py @@ -1,5 +1,4 @@ - -""" +r""" Represent hexagonal 2D mesh. The hexagons are pointy topped (as against flat topped). We use axial @@ -38,6 +37,7 @@ Coordinates a """ + import numpy as np from math import sqrt diff --git a/fidimag/atomistic/lib/clib.h b/fidimag/atomistic/lib/clib.h deleted file mode 100644 index 41fc895e..00000000 --- a/fidimag/atomistic/lib/clib.h +++ /dev/null @@ -1,122 +0,0 @@ -#ifndef __CLIB__ -#define __CLIB__ - -#include -#include -#include -//#include - -#include "fidimag_random.h" - -#define WIDE_PI 3.1415926535897932384626433832795L - -// ---------------------------------------------------------------------------- -/* 3 components for the cross product calculations */ - -inline double cross_x(double a0, double a1, double a2, double b0, double b1, - double b2) { - return a1 * b2 - a2 * b1; -} -inline double cross_y(double a0, double a1, double a2, double b0, double b1, - double b2) { - return a2 * b0 - a0 * b2; -} -inline double cross_z(double a0, double a1, double a2, double b0, double b1, - double b2) { - return a0 * b1 - a1 * b0; -} - -// ---------------------------------------------------------------------------- -// From exch.c - -void compute_exch_field(double *restrict spin, double *restrict field, double *restrict mu_s_inv, - double *restrict energy, double Jx, - double Jy, double Jz, int *restrict ngbs, int n, int n_ngbs); - -void compute_exch_field_spatial(double *restrict spin, double *restrict field, double *restrict mu_s_inv, - double *restrict energy, - double *restrict J, int *restrict ngbs, int n, int n_ngbs); - -double compute_exch_energy(double *restrict spin, double Jx, double Jy, double Jz, - int nx, int ny, int nz, int xperiodic, - int yperiodic); - -void compute_full_exch_field(double *restrict spin, double *restrict field, double *restrict mu_s_inv, - double *restrict energy, - double *restrict J, int *restrict ngbs, int n, int n_ngbs, - int n_shells, int *restrict n_ngbs_shell, int *restrict sum_ngbs_shell - ); - -// ----------------------------------------------------------------------------- -// From anis.c - -void compute_anis(double *restrict spin, double *restrict field, - double *restrict mu_s_inv, - double *restrict energy, double *restrict Ku, - double *restrict axis, int n); - -void compute_anis_cubic(double *restrict spin, double *restrict field, - double *restrict mu_s_inv, - double *restrict energy, - double *Kc, int n); - -// ---------------------------------------------------------------------------- -// From dmi.c - -void dmi_field_bulk(double *restrict spin, double *restrict field, - double *restrict mu_s_inv, - double *restrict energy, double *D, - int *restrict ngbs, int n, int n_ngbs); - -void dmi_field_interfacial_atomistic(double *spin, double *field, - double *mu_s_inv, - double *energy, double D, int *ngbs, int n, - int n_ngbs, int n_ngbs_dmi, double *DMI_vec); - -double dmi_energy(double *restrict spin, double D, int nx, int ny, int nz, int xperiodic, - int yperiodic); - -// ---------------------------------------------------------------------------- -// From demag_full.c - -void demag_full(double *restrict spin, double *restrict field, double *restrict energy, double *restrict coords, - double *restrict mu_s, double *restrict mu_s_scale, int n); - -// ---------------------------------------------------------------------------- -// From util.c - - -double skyrmion_number(double *spin, double *charge, int nx, int ny, int nz, - int *ngbs, int n_ngbs); - -double skyrmion_number_BergLuscher(double *spin, double *charge, int nx, int ny, - int nz, int *ngbs, int n_ngbs); - -void compute_guiding_center(double *spin, int nx, int ny, int nz, int nx_start, - int nx_stop, int ny_start, int ny_stop, - double *res); - -void compute_px_py_c(double *spin, int nx, int ny, int nz, double *px, - double *py); - -// ---------------------------------------------------------------------------- -// From sllg.c - -void normalise(double *m, int *pins, int n); - -void llg_rhs_dw_c(double *restrict m, double *restrict h, double *restrict dm, double *restrict T, double *restrict alpha, - double *restrict mu_s_inv, int *pins, double *restrict eta, int n, double gamma, - double dt); - -// ---------------------------------------------------------------------------- -// From mc.c - -void llg_s_rhs(double *restrict dm_dt, double *restrict spin, double *restrict h, double *restrict alpha, - double *restrict chi, double gamma, int n); - -void run_step_mc(mt19937_state *state, double *spin, double *new_spin, - int *ngbs, int *nngbs, int n_ngbs, double J, double J1, double D, - double D1, double *h, double Kc, int n, double T, - int hexagonal_mesh); - -#endif diff --git a/fidimag/atomistic/llg.py b/fidimag/atomistic/llg.py index 03b9454f..35b8e480 100644 --- a/fidimag/atomistic/llg.py +++ b/fidimag/atomistic/llg.py @@ -2,7 +2,7 @@ from __future__ import print_function # Use the common C library for the LLG equation in the atomistic case -import fidimag.extensions.common_clib as clib +import fidimag.extensions.c_clib as clib from .atomistic_driver import AtomisticDriver diff --git a/fidimag/atomistic/llg_stt.py b/fidimag/atomistic/llg_stt.py index 2d5c05c1..e495bd38 100644 --- a/fidimag/atomistic/llg_stt.py +++ b/fidimag/atomistic/llg_stt.py @@ -1,6 +1,5 @@ from __future__ import division - -import fidimag.extensions.common_clib as clib +import fidimag.extensions.c_clib as clib import numpy as np from .atomistic_driver import AtomisticDriver diff --git a/fidimag/atomistic/llg_stt_cpp.py b/fidimag/atomistic/llg_stt_cpp.py index 7f0b8996..30028aa0 100644 --- a/fidimag/atomistic/llg_stt_cpp.py +++ b/fidimag/atomistic/llg_stt_cpp.py @@ -2,7 +2,7 @@ import types -import fidimag.extensions.common_clib as clib +import fidimag.extensions.c_clib as clib import numpy as np from .atomistic_driver import AtomisticDriver diff --git a/fidimag/atomistic/monte_carlo.py b/fidimag/atomistic/monte_carlo.py index 9d5cc15d..5be33bcd 100644 --- a/fidimag/atomistic/monte_carlo.py +++ b/fidimag/atomistic/monte_carlo.py @@ -2,7 +2,7 @@ from __future__ import print_function import os import numpy as np -import fidimag.extensions.clib as clib +import fidimag.extensions.a_clib as clib from fidimag.common.fileio import DataSaver import fidimag.common.helper as helper import fidimag.common.constant as const diff --git a/fidimag/atomistic/sim.py b/fidimag/atomistic/sim.py index 59e60e59..3b6d755f 100644 --- a/fidimag/atomistic/sim.py +++ b/fidimag/atomistic/sim.py @@ -7,7 +7,7 @@ import fidimag.common.skyrmion_number import fidimag.common.helper as helper -import fidimag.extensions.clib as clib +import fidimag.extensions.a_clib as clib import numpy as np from fidimag.common.sim_base import SimBase diff --git a/fidimag/atomistic/sllg.py b/fidimag/atomistic/sllg.py index 323437f7..c3b4d905 100644 --- a/fidimag/atomistic/sllg.py +++ b/fidimag/atomistic/sllg.py @@ -1,6 +1,6 @@ from __future__ import division import numpy as np -import fidimag.extensions.clib as clib +import fidimag.extensions.a_clib as clib import fidimag.common.helper as helper import fidimag.common.constant as const diff --git a/fidimag/common/c_clib.pyx b/fidimag/common/c_clib.pyx new file mode 100644 index 00000000..1d0632c6 --- /dev/null +++ b/fidimag/common/c_clib.pyx @@ -0,0 +1,469 @@ +# distutils: language = c++ + +cimport numpy as np +import numpy as np + +# C++ modules +from libcpp cimport bool + +# ----------------------------------------------------------------------------- + +cdef extern from "c_clib.h": + + # From: llg.c + void llg_rhs(double * dm_dt, double * spin, + double *h, double *alpha, int *pins, + double gamma, int n, int do_precession, double default_c) + + + void llg_rhs_jtimes(double *jtn, double *m, double *h, + double *mp, double *hp, double *alpha, int *pins, + double gamma, int n, + int do_precession, double default_c) + + void compute_stt_field_c(double *spin, double *field, + double *jx, double *jy, double *jz, + double dx, double dy, double dz, int *ngbs, int n) + + # From: stt.c + void llg_stt_rhs(double *dm_dt, double *m, double *h, double *h_stt, + double *alpha,double beta, double u0, double gamma, int n) + + void llg_stt_cpp(double *dm_dt, double *m, double *h, double *p, + double *alpha, int *pins, double *a_J, + double beta, double gamma, int n) + + # ------------------------------------------------------------------------- + # From steepest_descent.c + + void sd_update_spin (double *spin, double *spin_last, double *magnetisation, + double *mxH, double *mxmxH, double *mxmxH_last, double tau, + int* pins, int n) + + void sd_compute_step (double *spin, double *spin_last, double *magnetisation, + double *field, + double *mxH, double *mxmxH, double *mxmxH_last, double tau, + int *pins, int n, int counter, double tmin, double tmax) + +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- + +def compute_llg_rhs(double [:] dm_dt, + double [:] spin, + double [:] field, + double [:] alpha, + int [:] pins, + gamma, n, do_precession, default_c): + llg_rhs(&dm_dt[0], &spin[0], &field[0], &alpha[0], &pins[0], + gamma, n, do_precession, default_c) + + +def compute_llg_jtimes(double [:] jtn, + double [:] m, + double [:] field, + double [:] mp, + double [:] field_p, + double [:] alpha, + int [:] pins, + gamma, n, do_precession, default_c): + llg_rhs_jtimes(&jtn[0], &m[0], &field[0], &mp[0], &field_p[0], + &alpha[0], &pins[0], gamma, n, do_precession, default_c) + +# ----------------------------------------------------------------------------- + +def compute_stt_field(double [:] spin, + double [:] field, + double [:] jx, + double [:] jy, + double [:] jz, + dx, dy, dz, + int [:, :] ngbs, + n + ): + compute_stt_field_c(&spin[0], &field[0], &jx[0], &jy[0],&jz[0], + dx, dy, dz, &ngbs[0, 0], n) + +def compute_llg_stt_rhs(double [:] dm_dt, + double [:] spin, + double [:] field, + double [:] field_stt, + double [:] alpha, + beta, u0, gamma, n): + llg_stt_rhs(&dm_dt[0], &spin[0], &field[0], &field_stt[0], + &alpha[0], beta, u0, gamma, n) + + + +def compute_llg_stt_cpp(double [:] dm_dt, + double [:] spin, + double [:] field, + double [:] p, + double [:] alpha, + int [:] pin, + double [:] a_J, + beta, gamma, n): + llg_stt_cpp(&dm_dt[0], &spin[0], &field[0], &p[0], + &alpha[0], &pin[0], &a_J[0], beta, gamma, n) + + +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- + + +def compute_sd_spin(double [:] spin, + double [:] spin_last, + double [:] magnetisation, + double [:] mxH, + double [:] mxmxH, + double [:] mxmxH_last, + double tau, + int [:] pins, + n): + + sd_update_spin(&spin[0], &spin_last[0], &magnetisation[0], &mxH[0], + &mxmxH[0], &mxmxH_last[0], tau, &pins[0], n + ) + +def compute_sd_step(double [:] spin, + double [:] spin_last, + double [:] magnetisation, + double [:] field, + double [:] mxH, + double [:] mxmxH, + double [:] mxmxH_last, + double tau, + int [:] pins, + n, counter, tmin, tmax): + + sd_compute_step(&spin[0], &spin_last[0], &magnetisation[0], + &field[0], &mxH[0], + &mxmxH[0], &mxmxH_last[0], tau, &pins[0], + n, counter, tmin, tmax + ) + +def normalise(a): + """ + normalise the given array a + """ + a.shape = (-1, 3) + b = np.sqrt(a[:, 0] ** 2 + a[:, 1] ** 2 + a[:, 2] ** 2) + ids = (b == 0) + b[ids] = 1.0 + a[:, 0] /= b + a[:, 1] /= b + a[:, 2] /= b + a.shape = (-1,) + +def init_scalar(value, mesh, *args): + + n = mesh.n + + mesh_v = np.zeros(n) + + if isinstance(value, (int, float)): + mesh_v[:] = value + elif hasattr(value, '__call__'): + for i in range(n): + mesh_v[i] = value(mesh.coordinates[i], *args) + + elif isinstance(value, np.ndarray): + if value.shape == mesh_v.shape: + mesh_v[:] = value[:] + else: + raise ValueError("Array size must match the mesh size") + + return mesh_v + +def init_vector(m0, mesh, dim=3, norm=False, *args): + n = mesh.n + spin = np.zeros((n, dim)) + if isinstance(m0, list) or isinstance(m0, tuple): + spin[:, :] = m0 + spin = np.reshape(spin, dim * n, order='C') + elif hasattr(m0, '__call__'): + v = m0(mesh.coordinates[0], *args) + if len(v) != dim: + raise Exception( + 'The length of the value in init_vector method must be {}.'.format(dim)) + for i in range(n): + spin[i, :] = m0(mesh.coordinates[i], *args) + spin = np.reshape(spin, dim * n, order='C') + elif isinstance(m0, np.ndarray): + if m0.shape == (dim, ): + spin[:] = m0 # broadcasting + else: + spin.shape = (-1) + spin[:] = m0 # overwriting the whole thing + spin.shape = (-1,) + if norm: + normalise(spin) + return spin + +def init_vector_func_fast(m0, mesh, double[:] field, norm=False, *args): + """ + + An unsafe method of setting the field. Depends on the setter code being + memory safe. + + m0 must be a Python function that takes the mesh and field as arguments. + Within that, the user must handle evaluating the function at different + coordinate points. It needs to be able to handle the spatial dependence + itself, and write the field valuse into the field array. This can be + written with Cython which will give much better performance. For example: + + from libc.math cimport sin + + def fast_sin_init(mesh, double[:] field, *params): + t, axis, Bmax, fc = params + for i in range(mesh.n): + field[3*i+0] = Bmax * axis[0] * sin(fc*t) + field[3*i+1] = Bmax * axis[1] * sin(fc*t) + field[3*i+2] = Bmax * axis[2] * sin(fc*t) + + """ + m0(mesh, field, *args) + if norm: + normalise(field) + return field + +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- +# CLASSES + +# ----------------------------------------------------------------------------- +# C++ definitions + +# See: +# https://cython.readthedocs.io/en/latest/src/userguide/wrapping_CPlusPlus.html +# +# - Methods that are redefined in inherited classes are only specified in the +# base class +# - These can be declared in a .pxd file: + +cdef extern from "c_energy.h": + + cdef cppclass Energy: + # except +: Without this declaration, C++ exceptions originating from + # the constructor will not be handled by Cython. + Energy() except + + + void compute_field(double t) + double compute_energy() + + # # Not using these variables from Cython: + # bool set_up + # int nx, ny, nz, n + # double dx, dy, dz + # double unit_length + # double *spin + # double *Ms + # double *Ms_inv + # double *field + # double *energy + # double *coordinates + # int *ngbs + int interaction_id + + cdef cppclass ExchangeEnergy(Energy): + ExchangeEnergy() except + + void setup(double * A, MicroSim * sim) + # double *A + + cdef cppclass AnisotropyEnergy(Energy): + AnisotropyEnergy() except + + void setup(double * Ku, double * axis, MicroSim * sim) + + cdef cppclass DMIEnergy(Energy): + ExchangeEnergy() except + + void setup(double * D, double * dmi_vector, int n_dmis, MicroSim * sim) + +# cdef extern from "c_energy.cpp": +# pass + +# ----------------------------------------------------------------------------- +# Python definitions + +# Cython initializes C++ class attributes of a cdef class using the nullary +# constructor. If the class you’re wrapping does not have a nullary +# constructor, you must store a pointer to the wrapped class and manually +# allocate and deallocate it. A convenient and safe place to do so is in the +# __cinit__ and __dealloc__ methods which are guaranteed to be called exactly +# once upon creation and deletion of the Python instance. + +cdef class PyEnergy: + cdef Energy *thisptr + def __cinit__(self): + # No need to allocate memory, we are not using the C++ base class + # self.thisptr = new Energy() + if type(self) != PyEnergy: + return + print("In Python A") + + # Necessary?: + def compute_field(self, t): + self.thisptr.compute_field(t) + + def compute_energy(self, time): + return self.thisptr.compute_energy() + + def add_interaction_to_sim(self, PyMicroSim sim): + sim.thisptr.add_interaction( self.thisptr) + + def get_interaction_id(self): + return self._thisptr.interaction_id + + +cdef class PyExchangeEnergy(PyEnergy): + cdef ExchangeEnergy *_thisptr + # Try cinit: + def __cinit__(self, double [:] A, PyMicroSim sim): + + self._thisptr = self.thisptr = new ExchangeEnergy() + self._thisptr.setup(&A[0], sim.thisptr) + + def setup(self, double [:] A, PyMicroSim sim): + self._thisptr.setup(&A[0], sim.thisptr) + + def __dealloc__(self): + del self._thisptr + + +cdef class PyAnisotropyEnergy(PyEnergy): + cdef AnisotropyEnergy *_thisptr + # Try cinit: + def __cinit__(self, double [:] Ku, double [:] axis, PyMicroSim sim): + + self._thisptr = self.thisptr = new AnisotropyEnergy() + self._thisptr.setup(&Ku[0], &axis[0], sim.thisptr) + + def setup(self, double [:] Ku, double [:] axis, PyMicroSim sim): + self._thisptr.setup(&Ku[0], &axis[0], sim.thisptr) + + def __dealloc__(self): + del self._thisptr + + +cdef class PyDMIEnergy(PyEnergy): + cdef DMIEnergy *_thisptr + # Try cinit: + def __cinit__(self, double [:] D, double [:] dmi_vector, n_dmis, + PyMicroSim sim): + + self._thisptr = self.thisptr = new DMIEnergy() + self._thisptr.setup(&D[0], &dmi_vector[0], n_dmis, sim.thisptr) + + def setup(self, double [:] D, double [:] dmi_vector, n_dmis, PyMicroSim sim): + self._thisptr.setup(&D[0], &dmi_vector[0], n_dmis, sim.thisptr) + + def __dealloc__(self): + del self._thisptr + + +# Simulation class ------------------------------------------------------------ + +cdef extern from "c_micro_sim.h": + + cdef cppclass MicroSim: + # except +: Without this declaration, C++ exceptions originating from + # the constructor will not be handled by Cython. + MicroSim() except + + + void setup(int nx, int ny, int nz, double dx, double dy, double dz, + double unit_length, double *coordinates, int *ngbs, + double *spin, double *Ms, double *Ms_inv, + double *energy, double *field, int *pins + ) + + void add_interaction(Energy * interaction) + + +cdef class PyMicroSim(object): + """ + Wrapper for the C++ MicroSim simulation class + """ + cdef MicroSim *thisptr + # Try cinit: + def __cinit__(self): + + self.thisptr = new MicroSim() + + def __dealloc__(self): + del self.thisptr + + def setup(self, nx, ny, nz, dx, dy, dz, unit_length, + double [:, :] coordinates, int [:, :] neighbours, + double [:] spin, double [:] Ms, double [:] Ms_inv, + double [:] energy, double [:] field, int [:] pins + ): + + return self.thisptr.setup(nx, ny, nz, dx, dy, dz, unit_length, + &coordinates[0, 0], &neighbours[0, 0], + &spin[0], &Ms[0], &Ms_inv[0], + &energy[0], &field[0], &pins[0] + ) + + def add_interaction(self, Interaction): + Interaction.add_interaction_to_sim(self) + +# Drivers --------------------------------------------------------------------- + +cdef extern from "m_driver.h": + + ctypedef struct LLG_params: + double gamma + double * alpha + double * field + int * pins + + + ctypedef enum IntegratorID: + RK4 + + cdef cppclass MicroLLGDriver: + # except +: Without this declaration, C++ exceptions originating from + # the constructor will not be handled by Cython. + MicroLLGDriver() except + + + void setup(MicroSim * sim, double * alpha, double gamma, double t, double dt) + void add_integrator(IntegratorID integrator_id) + + double t + double dt + + # LLG_params struct? + LLG_params llg_params + +cpdef enum: + INTEGRATOR_RK4 = 0 + +cdef class PyMicroLLGDriver(object): + """ + Wrapper for the C++ Micro LLG Driver simulation class + """ + cdef MicroLLGDriver *thisptr + # Try cinit: + def __cinit__(self): + + self.thisptr = new MicroLLGDriver() + + def __dealloc__(self): + del self.thisptr + + def setup(self, PyMicroSim sim, double [:] alpha, gamma, t, dt): + + return self.thisptr.setup( sim.thisptr, + &alpha[0], gamma, t, dt) + + def add_integrator(self, integrator_id): + return self.thisptr.add_integrator(integrator_id) + +# Test ------------------------------------------------------------------------ + +cdef extern from "m_driver_test.h": + + void test_integrator_RK4() + + +def C_test_integrator_RK4(): + + test_integrator_RK4() diff --git a/fidimag/common/chain_method_base.py b/fidimag/common/chain_method_base.py index 2463f935..b5af98e7 100644 --- a/fidimag/common/chain_method_base.py +++ b/fidimag/common/chain_method_base.py @@ -436,7 +436,9 @@ def initialise_integrator(self, self.n_dofs_image, mass=1, stepsize=1e-4) - self.integrator.set_options() + # Note: disabled as at the moment, set_options does not + # do anything for the Verlet integrator! + # self.integrator.set_options() # In Verlet algorithm we only use the total force G and not YxYxG: self._llg_evolve = False else: diff --git a/fidimag/common/chain_method_integrators.py b/fidimag/common/chain_method_integrators.py index 01c379bd..93010f3f 100644 --- a/fidimag/common/chain_method_integrators.py +++ b/fidimag/common/chain_method_integrators.py @@ -86,7 +86,7 @@ def run_until(self, t): return 0 def set_options(self, rtol=1e-8, atol=1e-8): - warnings.warn("Tolerances not available for VerletIntegrator") + raise NotImplementedError("Tolerances not available for VerletIntegrator") def _step(self, t, y, h, f): """ diff --git a/fidimag/common/dipolar/dipolar.h b/fidimag/common/dipolar/dipolar.h deleted file mode 100644 index 5d68e0f2..00000000 --- a/fidimag/common/dipolar/dipolar.h +++ /dev/null @@ -1,117 +0,0 @@ -#include -#include -#include -//#include - -#define WIDE_PI 3.1415926535897932384626433832795L - -inline double cross_x(double a0, double a1, double a2, double b0, double b1, double b2) { return a1*b2 - a2*b1; } -inline double cross_y(double a0, double a1, double a2, double b0, double b1, double b2) { return a2*b0 - a0*b2; } -inline double cross_z(double a0, double a1, double a2, double b0, double b1, double b2) { return a0*b1 - a1*b0; } - - -enum Type_Nij { - Tensor_xx, Tensor_yy, Tensor_zz, Tensor_xy, Tensor_xz, Tensor_yz -}; - -//========================================== -//used for demag - -typedef struct { - int nx; - int ny; - int nz; - double dx; - double dy; - double dz; - int lenx; - int leny; - int lenz; - - int total_length; - - //TODO: free tensors after obtaining NXX to save memory? - double *tensor_xx; - double *tensor_yy; - double *tensor_zz; - double *tensor_xy; - double *tensor_xz; - double *tensor_yz; - - //TODO: (1)using double, (2)using small size - fftw_complex *Nxx; //Nxx, Nxy .. are pure real now. - fftw_complex *Nyy; - fftw_complex *Nzz; - fftw_complex *Nxy; - fftw_complex *Nxz; - fftw_complex *Nyz; - - fftw_complex *Mx; - fftw_complex *My; - fftw_complex *Mz; - fftw_complex *Hx; - fftw_complex *Hy; - fftw_complex *Hz; - - double *mx; - double *my; - double *mz; - double *hx; - double *hy; - double *hz; - - //we need three plans - fftw_plan tensor_plan; - fftw_plan m_plan; - fftw_plan h_plan; - -} fft_demag_plan; - -fft_demag_plan *create_plan(void); -void finalize_plan(fft_demag_plan *restrict plan); -void init_plan(fft_demag_plan *plan, double dx, double dy, - double dz, int nx, int ny, int nz); -void compute_dipolar_tensors(fft_demag_plan *restrict plan); -void compute_demag_tensors(fft_demag_plan *restrict plan); -void create_fftw_plan(fft_demag_plan *restrict plan); - -void compute_demag_tensors_2dpbc(fft_demag_plan *restrict plan, double *restrict tensors, double pbc_2d_error, int sample_repeat_nx, int sample_repeat_ny, double dipolar_radius); -void fill_demag_tensors_c(fft_demag_plan *restrict plan, double *restrict tensors); - -void compute_fields(fft_demag_plan *restrict plan, double *restrict spin, double *mu_s, double *restrict field); -void exact_compute(fft_demag_plan *restrict plan, double *restrict spin, double *mu_s, double *restrict field); -double compute_demag_energy(fft_demag_plan *restrict plan, double *restrict spin, double *restrict mu_s, double *restrict field, double *restrict energy); - - - -//========================================================= -//========================================================= -//used for sode -typedef struct { - int nxyz; - - double dt; - double T; - double gamma; - double *mu_s; - double coeff; - double Q; - - double theta; - double theta1; - double theta2; - - double *dm1; - double *dm2; - double *eta; - -} ode_solver; - -void init_solver(ode_solver *s, double k_B, double theta, int nxyz, double dt, double gamma); -ode_solver *create_ode_plan(void); -void finalize_ode_plan(ode_solver *plan); -void run_step1(ode_solver *s, double *m, double *h, double *m_pred, double *T, - double *alpha, double *mu_s_inv, int *pins); -void run_step2(ode_solver *s, double *m_pred, double *h, double *m, double *T, - double *alpha, double *mu_s_inv, int *pins); - diff --git a/fidimag/common/dipolar/dipolar.pyx b/fidimag/common/dipolar/dipolar.pyx index b1d27d89..8fddc6fa 100644 --- a/fidimag/common/dipolar/dipolar.pyx +++ b/fidimag/common/dipolar/dipolar.pyx @@ -1,10 +1,12 @@ +# distutils: language = c++ + import numpy as np import numpy import ctypes cimport numpy as np np.import_array() -cdef extern from "dipolar.h": +cdef extern from "c_dipolar.h": # used for demag ctypedef struct fft_demag_plan: int nx, ny, nz @@ -184,7 +186,7 @@ cdef class FFTDemag(object): &field[0], &energy[0]) -cdef extern from "demagcoef.h": +cdef extern from "c_demagcoef.h": double CalculateSDA00(double x, double y, double z, double dx,double dy,double dz) double DemagNxxAsymptotic(double x, double y, double z, double dx,double dy,double dz) diff --git a/fidimag/common/gradient_descent.py b/fidimag/common/gradient_descent.py index c1ea8c2d..123eaeb1 100644 --- a/fidimag/common/gradient_descent.py +++ b/fidimag/common/gradient_descent.py @@ -1,6 +1,6 @@ from __future__ import division import numpy as np -import fidimag.extensions.clib as clib +import fidimag.extensions.a_clib as a_clib # import fidimag.common.constant as const from .minimiser_base import MinimiserBase diff --git a/fidimag/common/helper.py b/fidimag/common/helper.py index 42816440..20f44ee3 100644 --- a/fidimag/common/helper.py +++ b/fidimag/common/helper.py @@ -1,8 +1,8 @@ -from fidimag.extensions.common_clib import normalise -from fidimag.extensions.common_clib import init_scalar -from fidimag.extensions.common_clib import init_vector -from fidimag.extensions.common_clib import init_vector_func_fast -import fidimag.extensions.clib as clib +from fidimag.extensions.c_clib import normalise +from fidimag.extensions.c_clib import init_scalar +from fidimag.extensions.c_clib import init_vector +from fidimag.extensions.c_clib import init_vector_func_fast +import fidimag.extensions.a_clib as clib import numpy as np diff --git a/fidimag/common/lib/common_clib.h b/fidimag/common/lib/common_clib.h deleted file mode 100644 index d7a024ab..00000000 --- a/fidimag/common/lib/common_clib.h +++ /dev/null @@ -1,79 +0,0 @@ -#ifndef __CLIB__ -#define __CLIB__ - -#include -//#include -#define WIDE_PI 3.1415926535897932384626433832795L - -// ---------------------------------------------------------------------------- - -/* 3 components for the cross product calculations */ -inline double cross_x(double a0, double a1, double a2, double b0, double b1, - double b2) { - return a1 * b2 - a2 * b1; -} -inline double cross_y(double a0, double a1, double a2, double b0, double b1, - double b2) { - return a2 * b0 - a0 * b2; -} -inline double cross_z(double a0, double a1, double a2, double b0, double b1, - double b2) { - return a0 * b1 - a1 * b0; -} - -inline void normalise(double *m, int *pins, int n){ - int i, j, k; - double mm; - for (int id = 0; id < n; id++) { - i = 3*id; - j = i + 1; - k = j + 1; - - if (pins[id]>0) continue; - - mm = sqrt(m[i] * m[i] + m[j] * m[j] + m[k] * m[k]); - if(mm > 0) { - mm = 1 / mm; - m[i] *= mm; - m[j] *= mm; - m[k] *= mm; - } - } -} - -// ---------------------------------------------------------------------------- -// From: llg.c - -void llg_rhs(double *restrict dm_dt, double *restrict spin, double *restrict h, double *restrict alpha, int *restrict pins, - double gamma, int n, int do_precession, double default_c); - -void llg_rhs_jtimes(double *restrict jtn, double *restrict m, double *restrict h, double *restrict mp, double *restrict hp, - double *restrict alpha, int *restrict pins, double gamma, int n, - int do_precession, double default_c); - -// ---------------------------------------------------------------------------- -// From: stt.c - -void compute_stt_field_c(double *restrict spin, double *restrict field, double *restrict jx, double *restrict jy, - double *restrict jz, double dx, double dy, double dz, int *restrict ngbs, - int n); - -void llg_stt_rhs(double *restrict dm_dt, double *restrict m, double *restrict h, double *restrict h_stt, - double *restrict alpha, double beta, double u0, double gamma, int n); - -void llg_stt_cpp(double *restrict dm_dt, double *restrict m, double *restrict h, double *restrict p, double *restrict alpha, - int *restrict pins, double *restrict a_J, double beta, double gamma, int n); - -// ---------------------------------------------------------------------------- -// From steepest_descent.c - -void sd_update_spin (double *spin, double *spin_last, double *magnetisation, - double *mxH, double *mxmxH, double *mxmxH_last, double tau, - int* pins, int n); - -void sd_compute_step (double *spin, double *spin_last, double *magnetisation, double *field, - double *mxH, double *mxmxH, double *mxmxH_last, double tau, - int *pins, int n, int counter, double tmin, double tmax); - - -#endif diff --git a/fidimag/common/lib/common_clib.pyx b/fidimag/common/lib/common_clib.pyx deleted file mode 100644 index 07606162..00000000 --- a/fidimag/common/lib/common_clib.pyx +++ /dev/null @@ -1,224 +0,0 @@ -cimport numpy as np -import numpy as np - -# ----------------------------------------------------------------------------- - -cdef extern from "common_clib.h": - - # From: llg.c - void llg_rhs(double * dm_dt, double * spin, - double *h, double *alpha, int *pins, - double gamma, int n, int do_precession, double default_c) - - - void llg_rhs_jtimes(double *jtn, double *m, double *h, - double *mp, double *hp, double *alpha, int *pins, - double gamma, int n, - int do_precession, double default_c) - - void compute_stt_field_c(double *spin, double *field, - double *jx, double *jy, double *jz, - double dx, double dy, double dz, int *ngbs, int n) - - # From: stt.c - void llg_stt_rhs(double *dm_dt, double *m, double *h, double *h_stt, - double *alpha,double beta, double u0, double gamma, int n) - - void llg_stt_cpp(double *dm_dt, double *m, double *h, double *p, - double *alpha, int *pins, double *a_J, - double beta, double gamma, int n) - - # ------------------------------------------------------------------------- - # From steepest_descent.c - - void sd_update_spin (double *spin, double *spin_last, double *magnetisation, - double *mxH, double *mxmxH, double *mxmxH_last, double tau, - int* pins, int n) - - void sd_compute_step (double *spin, double *spin_last, double *magnetisation, - double *field, - double *mxH, double *mxmxH, double *mxmxH_last, double tau, - int *pins, int n, int counter, double tmin, double tmax) - -# ----------------------------------------------------------------------------- -# ----------------------------------------------------------------------------- - -def compute_llg_rhs(double [:] dm_dt, - double [:] spin, - double [:] field, - double [:] alpha, - int [:] pins, - gamma, n, do_precession, default_c): - llg_rhs(&dm_dt[0], &spin[0], &field[0], &alpha[0], &pins[0], - gamma, n, do_precession, default_c) - - -def compute_llg_jtimes(double [:] jtn, - double [:] m, - double [:] field, - double [:] mp, - double [:] field_p, - double [:] alpha, - int [:] pins, - gamma, n, do_precession, default_c): - llg_rhs_jtimes(&jtn[0], &m[0], &field[0], &mp[0], &field_p[0], - &alpha[0], &pins[0], gamma, n, do_precession, default_c) - -# ----------------------------------------------------------------------------- - -def compute_stt_field(double [:] spin, - double [:] field, - double [:] jx, - double [:] jy, - double [:] jz, - dx, dy, dz, - int [:, :] ngbs, - n - ): - compute_stt_field_c(&spin[0], &field[0], &jx[0], &jy[0],&jz[0], - dx, dy, dz, &ngbs[0, 0], n) - -def compute_llg_stt_rhs(double [:] dm_dt, - double [:] spin, - double [:] field, - double [:] field_stt, - double [:] alpha, - beta, u0, gamma, n): - llg_stt_rhs(&dm_dt[0], &spin[0], &field[0], &field_stt[0], - &alpha[0], beta, u0, gamma, n) - - - -def compute_llg_stt_cpp(double [:] dm_dt, - double [:] spin, - double [:] field, - double [:] p, - double [:] alpha, - int [:] pin, - double [:] a_J, - beta, gamma, n): - llg_stt_cpp(&dm_dt[0], &spin[0], &field[0], &p[0], - &alpha[0], &pin[0], &a_J[0], beta, gamma, n) - - -# ----------------------------------------------------------------------------- -# ----------------------------------------------------------------------------- - - -def compute_sd_spin(double [:] spin, - double [:] spin_last, - double [:] magnetisation, - double [:] mxH, - double [:] mxmxH, - double [:] mxmxH_last, - double tau, - int [:] pins, - n): - - sd_update_spin(&spin[0], &spin_last[0], &magnetisation[0], &mxH[0], - &mxmxH[0], &mxmxH_last[0], tau, &pins[0], n - ) - -def compute_sd_step(double [:] spin, - double [:] spin_last, - double [:] magnetisation, - double [:] field, - double [:] mxH, - double [:] mxmxH, - double [:] mxmxH_last, - double tau, - int [:] pins, - n, counter, tmin, tmax): - - sd_compute_step(&spin[0], &spin_last[0], &magnetisation[0], - &field[0], &mxH[0], - &mxmxH[0], &mxmxH_last[0], tau, &pins[0], - n, counter, tmin, tmax - ) - -def normalise(a): - """ - normalise the given array a - """ - a.shape = (-1, 3) - b = np.sqrt(a[:, 0] ** 2 + a[:, 1] ** 2 + a[:, 2] ** 2) - ids = (b == 0) - b[ids] = 1.0 - a[:, 0] /= b - a[:, 1] /= b - a[:, 2] /= b - a.shape = (-1,) - -def init_scalar(value, mesh, *args): - - n = mesh.n - - mesh_v = np.zeros(n) - - if isinstance(value, (int, float)): - mesh_v[:] = value - elif hasattr(value, '__call__'): - for i in range(n): - mesh_v[i] = value(mesh.coordinates[i], *args) - - elif isinstance(value, np.ndarray): - if value.shape == mesh_v.shape: - mesh_v[:] = value[:] - else: - raise ValueError("Array size must match the mesh size") - - return mesh_v - -def init_vector(m0, mesh, dim=3, norm=False, *args): - n = mesh.n - spin = np.zeros((n, dim)) - if isinstance(m0, list) or isinstance(m0, tuple): - spin[:, :] = m0 - spin = np.reshape(spin, dim * n, order='C') - elif hasattr(m0, '__call__'): - v = m0(mesh.coordinates[0], *args) - if len(v) != dim: - raise Exception( - 'The length of the value in init_vector method must be {}.'.format(dim)) - for i in range(n): - spin[i, :] = m0(mesh.coordinates[i], *args) - spin = np.reshape(spin, dim * n, order='C') - elif isinstance(m0, np.ndarray): - if m0.shape == (dim, ): - spin[:] = m0 # broadcasting - else: - spin.shape = (-1) - spin[:] = m0 # overwriting the whole thing - spin.shape = (-1,) - if norm: - normalise(spin) - return spin - -def init_vector_func_fast(m0, mesh, double[:] field, norm=False, *args): - """ - An unsafe method of setting the field. Depends on - the setter code being memory safe. - - m0 must be a Python function that takes the mesh - and field as arguments. Within that, the user - must handle evaluating the function at different - coordinate points. It needs to be able to handle - the spatial dependence itself, and write the - field valuse into the field array. This can - be written with Cython which will give much - better performance. For example: - - from libc.math cimport sin - - def fast_sin_init(mesh, double[:] field, *params): - t, axis, Bmax, fc = params - for i in range(mesh.n): - field[3*i+0] = Bmax * axis[0] * sin(fc*t) - field[3*i+1] = Bmax * axis[1] * sin(fc*t) - field[3*i+2] = Bmax * axis[2] * sin(fc*t) - - """ - m0(mesh, field, *args) - if norm: - normalise(field) - return field diff --git a/fidimag/common/neb_method/nebm_clib.pyx b/fidimag/common/neb_method/nebm_clib.pyx index 2bfeca0e..ccc3964e 100644 --- a/fidimag/common/neb_method/nebm_clib.pyx +++ b/fidimag/common/neb_method/nebm_clib.pyx @@ -1,4 +1,8 @@ -cdef extern from "nebm_lib.h": +# distutils: language = c++ + +cimport fidimag.common.vectormath as vmath + +cdef extern from "c_nebm_lib.h": void compute_tangents_C(double *tangents, double *y, double *energies, @@ -21,7 +25,7 @@ cdef extern from "nebm_lib.h": int n_images, int n_dofs_image) - void normalise(double * a, int n) + void normalise_images_C(double * y, int n_images, int n_dofs_image) void normalise_spins_C(double * y, int n_images, int n_dofs_image) double compute_distance_cartesian(double * A, double * B, int n_dofs_image, @@ -57,7 +61,7 @@ cdef extern from "nebm_lib.h": int n_images, int n_dofs_image) -cdef extern from "nebm_integrators.h": +cdef extern from "c_nebm_integrators.h": double step_Verlet_C(double * forces, double * forces_prev, double * velocities, @@ -110,7 +114,7 @@ def compute_effective_force(double [:] G, ) def normalise_clib(double [:] a, n): - normalise(&a[0], n) + vmath.normalise(&a[0], n) def project_images(double [:] vector, double [:] y, @@ -175,7 +179,7 @@ def compute_dYdt_nc(double [:] y, -cdef extern from "nebm_geodesic_lib.h": +cdef extern from "c_nebm_geodesic_lib.h": double compute_geodesic_GreatCircle(double * A, double * B, int n_dofs_image, int * material, @@ -228,7 +232,7 @@ def image_distances_GreatCircle(double [:] distances, ) -cdef extern from "nebm_spherical_lib.h": +cdef extern from "c_nebm_spherical_lib.h": void normalise_spherical(double * a, int n) void normalise_images_spherical_C(double * y, int n_images, int n_dofs_image) diff --git a/fidimag/common/neb_method/nebm_lib.h b/fidimag/common/neb_method/nebm_lib.h deleted file mode 100644 index 6ebd9699..00000000 --- a/fidimag/common/neb_method/nebm_lib.h +++ /dev/null @@ -1,72 +0,0 @@ -#include "math.h" -#define WIDE_PI 3.1415926535897932384626433832795L - - -void cross_product(double * output, double * A, double * B); - -double dot_product(double * A, double * B, int n); - -double compute_norm(double *restrict a, int n); - -void normalise(double *restrict a, int n); - -void compute_tangents_C(double *restrict ys, double *restrict energy, - double *restrict tangents, int image_num, int nodes - ); - -void compute_spring_force_C( - double *restrict spring_force, - double *restrict y, - double *restrict tangents, - double *restrict k, - int n_images, - int n_dofs_image, - double *restrict distances - ); - -void compute_effective_force_C(double *restrict G, - double *restrict tangents, - double *restrict gradientE, - double *restrict spring_force, - int *restrict climbing_image, - int n_images, - int n_dofs_image - ); - -void compute_image_distances(double *restrict distances, - double *restrict path_distances, - double *restrict y, - int n_images, - int n_dofs_image, - double (* compute_distance)(double *, - double *, - int, - int *, - int), - int *restrict material, - int n_dofs_image_material - ); - -void normalise_images_C(double *restrict y, int n_images, int n_dofs_image); - -void normalise_spins_C(double *restrict y, int n_images, int n_dofs_image); - -void project_images_C(double *restrict vector, double *restrict y, - int n_images, int n_dofs_image - ); - -void project_vector_C(double *restrict vector, double *restrict y, - int n_dofs_image - ); - -double compute_distance_cartesian(double *restrict A, double *restrict B, int n_dofs_image, - int *restrict material, int n_dofs_image_material - ); - -void compute_dYdt_C(double *restrict y, double *restrict G, double *restrict dYdt, - int *restrict pins, - int n_images, int n_dofs_image); - -void compute_dYdt_nc_C(double *restrict y, double *restrict G, double *restrict dYdt, - int *restrict pins, - int n_images, int n_dofs_image); diff --git a/fidimag/common/neb_method/nebm_spherical_lib.h b/fidimag/common/neb_method/nebm_spherical_lib.h deleted file mode 100644 index 3e46ffce..00000000 --- a/fidimag/common/neb_method/nebm_spherical_lib.h +++ /dev/null @@ -1,11 +0,0 @@ -#include "math.h" -#define WIDE_PI 3.1415926535897932384626433832795L - -void normalise_spherical(double *restrict a, int n); - -void normalise_images_spherical_C(double *restrict y, int n_images, - int n_dofs_image); - -double compute_distance_spherical(double *restrict A, double * B, int n, - int *restrict material, int n_dofs_image_material - ); diff --git a/fidimag/common/nebm_spherical.py b/fidimag/common/nebm_spherical.py index 9fd7e300..f20e0ef3 100644 --- a/fidimag/common/nebm_spherical.py +++ b/fidimag/common/nebm_spherical.py @@ -13,7 +13,7 @@ class NEBM_Spherical(ChainMethodBase): - """ + r""" ARGUMENTS ----------------------------------------------------------------- diff --git a/fidimag/common/plot.py b/fidimag/common/plot.py index 4d2154e0..322ec92a 100644 --- a/fidimag/common/plot.py +++ b/fidimag/common/plot.py @@ -228,12 +228,12 @@ def plot_micro(sim, component='all', filename=None, figsize=(10, 5), if cbar is True: if component is 'angle': - # Some special handling to print \pi + # Some special handling to print pi # rather than the numbers! cbar = ax.cax.colorbar(im, ticks=np.linspace(0, 2*np.pi, ncbarticks)) - cbarlabels = ['${:.1f} \pi$'.format(x/(np.pi)) + cbarlabels = [r'${:.1f} \pi$'.format(x/(np.pi)) if x != 0.0 else '0.0' for x in np.linspace(0, 2*np.pi, ncbarticks)] cbar.ax.set_yticklabels(cbarlabels) @@ -461,12 +461,12 @@ def plot_atom_cub(sim, component='all', filename=None, figsize=(10, 5), if cbar is True: if component is 'angle': - # Some special handling to print \pi + # Some special handling to print pi # rather than the numbers! cbar = ax.cax.colorbar(im, ticks=np.linspace(0, 2*np.pi, ncbarticks)) - cbarlabels = ['${:.1f} \pi$'.format(x/(np.pi)) + cbarlabels = [r'${:.1f} \pi$'.format(x/(np.pi)) if x != 0.0 else '0.0' for x in np.linspace(0, 2*np.pi, ncbarticks)] cbar.ax.set_yticklabels(cbarlabels) @@ -714,7 +714,7 @@ def plot_atom_hex(sim, component='all', filename=None, figsize=(10, 5), cbar = ax.cax.colorbar(im, ticks=np.linspace(0, 2*np.pi, ncbarticks)) - cbarlabels = ['${:.1f} \pi$'.format(x/(np.pi)) + cbarlabels = [r'${:.1f} \pi$'.format(x/(np.pi)) if x != 0.0 else '0.0' for x in np.linspace(0, 2*np.pi, ncbarticks)] cbar.ax.set_yticklabels(cbarlabels) diff --git a/fidimag/common/plt_helper.py b/fidimag/common/plt_helper.py index 896246a7..c33e0cdf 100644 --- a/fidimag/common/plt_helper.py +++ b/fidimag/common/plt_helper.py @@ -1,11 +1,8 @@ from __future__ import division -# import matplotlib as mpl -# mpl.use("Agg") -import fidimag.extensions.clib as clib +import fidimag.extensions.a_clib as clib import numpy as np import matplotlib.pyplot as plt -# from mpl_toolkits.mplot3d import Axes3D from matplotlib.colors import colorConverter from matplotlib.collections import PolyCollection diff --git a/fidimag/common/simple_minimiser.py b/fidimag/common/simple_minimiser.py index 995831bf..9607eb89 100644 --- a/fidimag/common/simple_minimiser.py +++ b/fidimag/common/simple_minimiser.py @@ -1,6 +1,6 @@ from __future__ import division import numpy as np -import fidimag.extensions.clib as clib +import fidimag.extensions.a_clib as clib # import fidimag.common.constant as const from .minimiser_base import MinimiserBase diff --git a/fidimag/common/skyrmion_number.py b/fidimag/common/skyrmion_number.py index 992d5a48..fec66a61 100644 --- a/fidimag/common/skyrmion_number.py +++ b/fidimag/common/skyrmion_number.py @@ -1,4 +1,4 @@ -import fidimag.extensions.clib +import fidimag.extensions.a_clib import fidimag.extensions.micro_clib import numpy as np @@ -69,7 +69,7 @@ def skyrmion_number_slice(sim, at=None, zIndex=None): spinSlice, sim._skx_number, sim.mesh.nx, sim.mesh.ny, sim.mesh.nz, sim.mesh.neighbours) else: - return fidimag.extensions.clib.compute_skyrmion_number(\ + return fidimag.extensions.a_clib.compute_skyrmion_number(\ spinSlice, sim._skx_number, sim.mesh.nx, sim.mesh.ny, sim.mesh.nz, sim.mesh.neighbours) @@ -104,7 +104,7 @@ def skyrmion_number_lee(sim): sim.mesh.ny, sim.mesh.nz, sim.mesh.neighbours) else: - skyrmionNumbers[zI] = fidimag.extensions.clib.compute_skyrmion_number(\ + skyrmionNumbers[zI] = fidimag.extensions.a_clib.compute_skyrmion_number(\ spinSlice, sim._skx_number, sim.mesh.nx, sim.mesh.ny, sim.mesh.nz, sim.mesh.neighbours, 6) diff --git a/fidimag/common/steepest_descent.py b/fidimag/common/steepest_descent.py index ddbe0004..e1c0c16e 100644 --- a/fidimag/common/steepest_descent.py +++ b/fidimag/common/steepest_descent.py @@ -1,8 +1,8 @@ from __future__ import division import numpy as np -import fidimag.extensions.common_clib as clib +import fidimag.extensions.c_clib as clib # Change int he future to common clib: -import fidimag.extensions.clib as atom_clib +import fidimag.extensions.a_clib as atom_clib import sys from .minimiser_base import MinimiserBase diff --git a/fidimag/common/sundials/cvode.pyx b/fidimag/common/sundials/cvode.pyx index 92b76abc..2db74bfd 100644 --- a/fidimag/common/sundials/cvode.pyx +++ b/fidimag/common/sundials/cvode.pyx @@ -1,3 +1,5 @@ +# distutils: language = c++ + import numpy as np import os cimport numpy as np # import special compile-time information about numpy @@ -5,6 +7,7 @@ cimport openmp np.import_array() # don't remove or you'll segfault from libc.string cimport memcpy import sys +from fidimag.common cimport vectormath cdef extern from "sundials/sundials_types.h": @@ -14,7 +17,7 @@ cdef extern from "sundials/sundials_types.h": cdef extern from "sundials/sundials_nvector.h": cdef struct _generic_N_Vector: void *content - + ctypedef _generic_N_Vector *N_Vector N_Vector N_VNew_Serial(long int vec_length) N_Vector N_VNew_OpenMP(long int vec_length, int num_threads) @@ -168,16 +171,10 @@ cdef extern from "sundials/sundials_iterative.h": int PREC_LEFT int PREC_RIGHT int PREC_BOTH - + int MODIFIED_GS int CLASSICAL_GS - - -cdef extern from "../../atomistic/lib/clib.h": - void normalise(double * m, int nxyz) - - cdef struct cv_userdata: void * rhs_fun void * y diff --git a/fidimag/common/vectormath.c b/fidimag/common/vectormath.c new file mode 100644 index 00000000..06f2230d --- /dev/null +++ b/fidimag/common/vectormath.c @@ -0,0 +1 @@ +#error Do not use this file, it is the result of a failed Cython compilation. diff --git a/fidimag/common/vectormath.pxd b/fidimag/common/vectormath.pxd new file mode 100644 index 00000000..ec9fbd73 --- /dev/null +++ b/fidimag/common/vectormath.pxd @@ -0,0 +1,3 @@ +cdef extern from "c_vectormath.h": + void normalise(double *m, int *pins, int n) + void normalise(double * a, int n) \ No newline at end of file diff --git a/fidimag/micro/exchange.py b/fidimag/micro/exchange.py index 5f6706b7..57864afb 100644 --- a/fidimag/micro/exchange.py +++ b/fidimag/micro/exchange.py @@ -1,6 +1,5 @@ import fidimag.extensions.micro_clib as micro_clib from .energy import Energy -#from constant import mu_0 class UniformExchange(Energy): @@ -9,12 +8,12 @@ class UniformExchange(Energy): UniformExchange(A, name='UniformExchange') Compute the exchange field in micromagnetics. - + Inputs: A: float - A is the exchange stiffness constant measured in + A is the exchange stiffness constant measured in Joules / Meter (J / M) - + """ def __init__(self, A, name='UniformExchange'): diff --git a/fidimag/micro/lib/anis.c b/fidimag/micro/lib/anis.c deleted file mode 100644 index f4966e7e..00000000 --- a/fidimag/micro/lib/anis.c +++ /dev/null @@ -1,32 +0,0 @@ -#include "micro_clib.h" - - -void compute_uniaxial_anis(double *restrict m, double *restrict field, double *restrict energy, double *restrict Ms_inv, - double *restrict Ku, double *restrict axis, int nx, int ny, int nz) { - - int n = nx * ny * nz; - - #pragma omp parallel for - for (int i = 0; i < n; i++) { - int j = 3 * i; - - if (Ms_inv[i] == 0.0){ - field[j] = 0; - field[j + 1] = 0; - field[j + 2] = 0; - energy[i] = 0; - continue; - } - - double m_u = m[j] * axis[j] + m[j + 1] * axis[j + 1] + m[j + 2] * axis[j + 2]; - - field[j] = 2 * Ku[i] * m_u * Ms_inv[i] * MU0_INV * axis[j]; - field[j + 1] = 2 * Ku[i] * m_u * Ms_inv[i] * MU0_INV * axis[j + 1]; - field[j + 2] = 2 * Ku[i] * m_u * Ms_inv[i] * MU0_INV * axis[j + 2]; - - energy[i] = Ku[i] * (1 - m_u * m_u); - - } - -} - diff --git a/fidimag/micro/lib/baryakhtar/baryakhtar_clib.pyx b/fidimag/micro/lib/baryakhtar_clib.pyx similarity index 97% rename from fidimag/micro/lib/baryakhtar/baryakhtar_clib.pyx rename to fidimag/micro/lib/baryakhtar_clib.pyx index 803b91b1..e06aeaf6 100644 --- a/fidimag/micro/lib/baryakhtar/baryakhtar_clib.pyx +++ b/fidimag/micro/lib/baryakhtar_clib.pyx @@ -1,8 +1,9 @@ +# distutils: language = c++ import numpy cimport numpy as np np.import_array() -cdef extern from "baryakhtar_clib.h": +cdef extern from "m_baryakhtar_clib.h": void compute_laplace_m(double *m, double *field, double *Ms, double dx, double dy, double dz, int nx, int ny, int nz) diff --git a/fidimag/micro/lib/micro_clib.h b/fidimag/micro/lib/micro_clib.h deleted file mode 100644 index 17cac732..00000000 --- a/fidimag/micro/lib/micro_clib.h +++ /dev/null @@ -1,32 +0,0 @@ -#include - -#include - -#define WIDE_PI 3.1415926535897932384626433832795L -#define MU0 1.25663706143591728850e-6 -#define MU0_INV 795774.71545947669074 - -inline double cross_x(double a0, double a1, double a2, - double b0, double b1, double b2) { return a1*b2 - a2*b1; } -inline double cross_y(double a0, double a1, double a2, - double b0, double b1, double b2) { return a2*b0 - a0*b2; } -inline double cross_z(double a0, double a1, double a2, - double b0, double b1, double b2) { return a0*b1 - a1*b0; } - -void compute_exch_field_micro(double *restrict m, double *restrict field, double *restrict energy, double *restrict Ms_inv, - double A, double dx, double dy, double dz, int n, int *ngbs); - -void dmi_field(double *restrict m, double *restrict field, - double *restrict energy, double *restrict Ms_inv, - double *restrict D, int n_DMIs, - double *dmi_vector, - double dx, double dy, double dz, int n, int *ngbs); - -void compute_exch_field_rkky_micro(double *m, double *field, double *energy, double *Ms_inv, - double sigma, int nx, double ny, double nz, int z_bottom, int z_top); - -void compute_uniaxial_anis(double *restrict m, double *restrict field, double *restrict energy, double *restrict Ms_inv, - double *restrict Ku, double *restrict axis, int nx, int ny, int nz); - -double skyrmion_number(double *restrict spin, double *restrict charge, - int nx, int ny, int nz, int *restrict ngbs); diff --git a/fidimag/micro/lib/micro_clib.pyx b/fidimag/micro/lib/micro_clib.pyx index 2cd3036a..72d6883b 100644 --- a/fidimag/micro/lib/micro_clib.pyx +++ b/fidimag/micro/lib/micro_clib.pyx @@ -1,7 +1,8 @@ +# distutils: language = c++ import numpy as np cimport numpy as np -cdef extern from "micro_clib.h": +cdef extern from "m_clib.h": void compute_exch_field_micro(double *m, double *field, double *energy, double *Ms_inv, double A, double dx, double dy, double dz, diff --git a/fidimag/micro/llg.py b/fidimag/micro/llg.py index 4e7f0c9d..15ad7e08 100644 --- a/fidimag/micro/llg.py +++ b/fidimag/micro/llg.py @@ -3,7 +3,7 @@ # We use the common/lib/llg.c file to calculate the LLG equation for # the micromagnetic case -import fidimag.extensions.common_clib as clib +import fidimag.extensions.c_clib as clib from .micro_driver import MicroDriver diff --git a/fidimag/micro/llg_stt.py b/fidimag/micro/llg_stt.py index 41fe7fba..021ec54f 100644 --- a/fidimag/micro/llg_stt.py +++ b/fidimag/micro/llg_stt.py @@ -1,6 +1,6 @@ from __future__ import division -import fidimag.extensions.common_clib as clib +import fidimag.extensions.c_clib as clib import numpy as np from .micro_driver import MicroDriver diff --git a/fidimag/micro/llg_stt_cpp.py b/fidimag/micro/llg_stt_cpp.py index 9b6c93cf..d272cba4 100644 --- a/fidimag/micro/llg_stt_cpp.py +++ b/fidimag/micro/llg_stt_cpp.py @@ -2,7 +2,7 @@ from .micro_driver import MicroDriver -import fidimag.extensions.common_clib as clib +import fidimag.extensions.c_clib as clib import numpy as np import fidimag.common.helper as helper diff --git a/fidimag/user/example/example.pyx b/fidimag/user/example/example.pyx index 29eeb3c7..431ab37f 100644 --- a/fidimag/user/example/example.pyx +++ b/fidimag/user/example/example.pyx @@ -1,3 +1,5 @@ +# distutils: language = c++ + from libc.math cimport cos, sin def fast_sin_init(mesh, double[:] field, *params): diff --git a/native/include/a_clib.h b/native/include/a_clib.h new file mode 100644 index 00000000..6a70f923 --- /dev/null +++ b/native/include/a_clib.h @@ -0,0 +1,101 @@ +#pragma once + +#include +#include +#include +#include +#include "a_random.h" +#include "c_vectormath.h" + + + +// ---------------------------------------------------------------------------- +// From exch.c + +void compute_exch_field(double *spin, double *field, double *mu_s_inv, + double *energy, double Jx, + double Jy, double Jz, int *ngbs, int n, int n_ngbs); + +void compute_exch_field_spatial(double * spin, double * field, double * mu_s_inv, + double * energy, + double * J, int * ngbs, int n, int n_ngbs); + +double compute_exch_energy(double * spin, double Jx, double Jy, double Jz, + int nx, int ny, int nz, int xperiodic, + int yperiodic); + +void compute_full_exch_field(double * spin, double * field, double * mu_s_inv, + double * energy, + double * J, int * ngbs, int n, int n_ngbs, + int n_shells, int * n_ngbs_shell, int * sum_ngbs_shell + ); + +// ----------------------------------------------------------------------------- +// From anis.c + +void compute_anis(double * spin, double * field, + double * mu_s_inv, + double * energy, double * Ku, + double * axis, int n); + +void compute_anis_cubic(double * spin, double * field, + double * mu_s_inv, + double * energy, + double *Kc, int n); + +// ---------------------------------------------------------------------------- +// From dmi.c + +void dmi_field_bulk(double * spin, double * field, + double * mu_s_inv, + double * energy, double *D, + int * ngbs, int n, int n_ngbs); + +void dmi_field_interfacial_atomistic(double *spin, double *field, + double *mu_s_inv, + double *energy, double D, int *ngbs, int n, + int n_ngbs, int n_ngbs_dmi, double *DMI_vec); + +double dmi_energy(double * spin, double D, int nx, int ny, int nz, int xperiodic, + int yperiodic); + +// ---------------------------------------------------------------------------- +// From demag_full.c + +void demag_full(double * spin, double * field, double * energy, double * coords, + double * mu_s, double * mu_s_scale, int n); + +// ---------------------------------------------------------------------------- +// From util.c + + +double skyrmion_number(double *spin, double *charge, int nx, int ny, int nz, + int *ngbs, int n_ngbs); + +double skyrmion_number_BergLuscher(double *spin, double *charge, int nx, int ny, + int nz, int *ngbs, int n_ngbs); + +void compute_guiding_center(double *spin, int nx, int ny, int nz, int nx_start, + int nx_stop, int ny_start, int ny_stop, + double *res); + +void compute_px_py_c(double *spin, int nx, int ny, int nz, double *px, + double *py); + +// ---------------------------------------------------------------------------- +// From sllg.c + +void llg_rhs_dw_c(double * m, double * h, double * dm, double * T, double * alpha, + double * mu_s_inv, int *pins, double * eta, int n, double gamma, + double dt); + +// ---------------------------------------------------------------------------- +// From mc.c + +void llg_s_rhs(double * dm_dt, double * spin, double * h, double * alpha, + double * chi, double gamma, int n); + +void run_step_mc(mt19937_state *state, double *spin, double *new_spin, + int *ngbs, int *nngbs, int n_ngbs, double J, double J1, double D, + double D1, double *h, double Kc, int n, double T, + int hexagonal_mesh); diff --git a/fidimag/atomistic/lib/fidimag_random.h b/native/include/a_random.h similarity index 83% rename from fidimag/atomistic/lib/fidimag_random.h rename to native/include/a_random.h index 64697cfa..ce038452 100644 --- a/fidimag/atomistic/lib/fidimag_random.h +++ b/native/include/a_random.h @@ -1,11 +1,4 @@ -#ifndef __FIDIMAG_RANDOM__ -#define __FIDIMAG_RANDOM__ - - - -//#include - -#define WIDE_PI 3.1415926535897932384626433832795L +#pragma once //================================================= //random number, mt19937 @@ -14,7 +7,7 @@ typedef struct { unsigned int matrix[2];// = { 0, 0x9908b0dfU}; int index_t; int seed; - + } mt19937_state; #define MT19973_RAND_MAX 4294967295u @@ -29,5 +22,3 @@ int rand_int_n(mt19937_state *state, int n); void gauss_random_vector(mt19937_state *state, double *x, int n); void uniform_random_sphere(mt19937_state *state, double *spin, int n); - -#endif diff --git a/native/include/c_clib.h b/native/include/c_clib.h new file mode 100644 index 00000000..3d2d9c14 --- /dev/null +++ b/native/include/c_clib.h @@ -0,0 +1,38 @@ +#pragma once +#include +#include + + +// ---------------------------------------------------------------------------- +// From: llg.c + +void llg_rhs(double * dm_dt, double * spin, double * h, double * alpha, int * pins, + double gamma, int n, int do_precession, double default_c); + +void llg_rhs_jtimes(double * jtn, double * m, double * h, double * mp, double * hp, + double * alpha, int * pins, double gamma, int n, + int do_precession, double default_c); + +// ---------------------------------------------------------------------------- +// From: stt.c + +void compute_stt_field_c(double * spin, double * field, double * jx, double * jy, + double * jz, double dx, double dy, double dz, int * ngbs, + int n); + +void llg_stt_rhs(double * dm_dt, double * m, double * h, double * h_stt, + double * alpha, double beta, double u0, double gamma, int n); + +void llg_stt_cpp(double * dm_dt, double * m, double * h, double * p, double * alpha, + int * pins, double * a_J, double beta, double gamma, int n); + +// ---------------------------------------------------------------------------- +// From steepest_descent.c + +void sd_update_spin (double *spin, double *spin_last, double *magnetisation, + double *mxH, double *mxmxH, double *mxmxH_last, double tau, + int* pins, int n); + +void sd_compute_step (double *spin, double *spin_last, double *magnetisation, double *field, + double *mxH, double *mxmxH, double *mxmxH_last, double tau, + int *pins, int n, int counter, double tmin, double tmax); diff --git a/native/include/c_constants.h b/native/include/c_constants.h new file mode 100644 index 00000000..a0f68537 --- /dev/null +++ b/native/include/c_constants.h @@ -0,0 +1,11 @@ +#pragma once + +// GNU compilers on Linux should already have this defined +#ifndef M_PIl +#define M_PIl 3.141592653589793238462643383279 +#endif + +#define MU_0 (4*M_PIl*1e-7) +#define MU_0_INV (1/MU_0) + + diff --git a/fidimag/common/dipolar/demagcoef.h b/native/include/c_demagcoef.h similarity index 84% rename from fidimag/common/dipolar/demagcoef.h rename to native/include/c_demagcoef.h index 94e34483..e56e4eaf 100644 --- a/fidimag/common/dipolar/demagcoef.h +++ b/native/include/c_demagcoef.h @@ -1,9 +1,10 @@ +#pragma once /* This file demag_oommf.h is taken from the OOMMF project (oommf/app/oxs/ext/demagcoef.h downloaded from http://math.nist.gov/oommf/dist/oommf12a5rc_20120928.tar.gz) -with slightly modifications (change OC_REALWIDE to double) +with slightly modifications (change OC_REALWIDE to double) and distributed under this license shown below. */ -/* License +/* License OOMMF - Object Oriented MicroMagnetic Framework @@ -130,49 +131,10 @@ double DemagNyzAsymptotic(double x, double y, double z, double dx,double dy,double dz); - //////////////////////////////////////////////////////////////////////////// // Routines to do accurate summation -static int AS_Compare(const void* px,const void* py) -{ - // Comparison based on absolute values - double x=fabs(*((const double *)px)); - double y=fabs(*((const double *)py)); - if(xy) return -1; - return 0; -} - - -static double -AccurateSum(int n,double *arr) -{ - // Order by decreasing magnitude - qsort(arr,n,sizeof(double),AS_Compare); - - // Add up using doubly compensated summation. If necessary, mark - // variables these "volatile" to protect against problems arising - // from extra precision. Also, don't expect the compiler to respect - // order with respect to parentheses at high levels of optimization, - // i.e., write "u=x; u-=(y-corr)" as opposed to "u=x-(y-corr)". - - double sum,corr,y,u,t,v,z,x,tmp; - - sum=arr[0]; corr=0; - for(int i=1;i +#include +#include +#include "c_vectormath.h" + +enum Type_Nij { + Tensor_xx, Tensor_yy, Tensor_zz, Tensor_xy, Tensor_xz, Tensor_yz +}; + +//========================================== +//used for demag + +typedef struct { + int nx; + int ny; + int nz; + double dx; + double dy; + double dz; + int lenx; + int leny; + int lenz; + + int total_length; + + //TODO: free tensors after obtaining NXX to save memory? + double *tensor_xx; + double *tensor_yy; + double *tensor_zz; + double *tensor_xy; + double *tensor_xz; + double *tensor_yz; + + //TODO: (1)using double, (2)using small size + std::complex *Nxx; //Nxx, Nxy .. are pure real now. + std::complex *Nyy; + std::complex *Nzz; + std::complex *Nxy; + std::complex *Nxz; + std::complex *Nyz; + + std::complex *Mx; + std::complex *My; + std::complex *Mz; + std::complex*Hx; + std::complex *Hy; + std::complex *Hz; + + double *mx; + double *my; + double *mz; + double *hx; + double *hy; + double *hz; + + //we need three plans + fftw_plan tensor_plan; + fftw_plan m_plan; + fftw_plan h_plan; + +} fft_demag_plan; + +fft_demag_plan *create_plan(void); +void finalize_plan(fft_demag_plan * plan); +void init_plan(fft_demag_plan *plan, double dx, double dy, + double dz, int nx, int ny, int nz); +void compute_dipolar_tensors(fft_demag_plan * plan); +void compute_demag_tensors(fft_demag_plan * plan); +void create_fftw_plan(fft_demag_plan * plan); + +void compute_demag_tensors_2dpbc(fft_demag_plan * plan, double * tensors, double pbc_2d_error, int sample_repeat_nx, int sample_repeat_ny, double dipolar_radius); +void fill_demag_tensors_c(fft_demag_plan * plan, double * tensors); + +void compute_fields(fft_demag_plan * plan, double * spin, double *mu_s, double * field); +void exact_compute(fft_demag_plan * plan, double * spin, double *mu_s, double * field); +double compute_demag_energy(fft_demag_plan * plan, double * spin, double * mu_s, double * field, double * energy); + + + +//========================================================= +//========================================================= +//used for sode +typedef struct { + int nxyz; + + double dt; + double T; + double gamma; + double *mu_s; + double coeff; + double Q; + + double theta; + double theta1; + double theta2; + + double *dm1; + double *dm2; + double *eta; + +} ode_solver; + +void init_solver(ode_solver *s, double k_B, double theta, int nxyz, double dt, double gamma); +ode_solver *create_ode_plan(void); +void finalize_ode_plan(ode_solver *plan); +void run_step1(ode_solver *s, double *m, double *h, double *m_pred, double *T, + double *alpha, double *mu_s_inv, int *pins); +void run_step2(ode_solver *s, double *m_pred, double *h, double *m, double *T, + double *alpha, double *mu_s_inv, int *pins); diff --git a/native/include/c_energy.h b/native/include/c_energy.h new file mode 100644 index 00000000..474ca676 --- /dev/null +++ b/native/include/c_energy.h @@ -0,0 +1,80 @@ +#pragma once +#include +#include +#include "c_micro_sim.h" + +class Energy { +public: + Energy() {}; + virtual ~Energy() {std::cout << "Killing A\n";}; + int interaction_id; + bool set_up; + int nx, ny, nz, n; + double dx, dy, dz; + double unit_length; + double *spin; + double *Ms; + double *Ms_inv; + double *field; + double *energy; + double *coordinates; + int *ngbs; + double compute_energy(); + // Will get the parameters from a simulation class + void _setup(MicroSim * sim); + virtual void compute_field(double t) {}; +}; + + +class ExchangeEnergy : public Energy { +public: + ExchangeEnergy() { + std::cout << "Instatiating ExchangeEnergy class; at " << this << "\n"; + this->interaction_id = 1; + }; + ~ExchangeEnergy() {std::cout << "Killing Exchange\n";}; + double *A; + void setup(double *A, MicroSim * sim) { + _setup(sim); + this->A = A; + } + void compute_field(double t); +}; + + +class AnisotropyEnergy : public Energy { +public: + AnisotropyEnergy() { + std::cout << "Instatiating AnisotropyEnergy class; at " << this << "\n"; + this->interaction_id = UNIAXIAL_ANISOTROPY; + }; + ~AnisotropyEnergy() {std::cout << "Killing Anisotropy\n";}; + double *Ku; + double *axis; + void setup(double *Ku, double *axis, MicroSim * sim) { + _setup(sim); + this->Ku = Ku; + this->axis = axis; + } + void compute_field(double t); +}; + + +class DMIEnergy : public Energy { +public: + DMIEnergy() { + std::cout << "Instatiating DMIEnergy class; at " << this << "\n"; + this->interaction_id = DMI; + }; + ~DMIEnergy() {std::cout << "Killing DMI\n";}; + double *D; + double *dmi_vector; + int n_dmis; + void setup(double * D, double *dmi_vector, int n_dmis, MicroSim * sim) { + _setup(sim); + this->D = D; + this->dmi_vector = dmi_vector; + this->n_dmis = n_dmis; + } + void compute_field(double t); +}; diff --git a/native/include/c_micro_sim.h b/native/include/c_micro_sim.h new file mode 100644 index 00000000..44fbb1d1 --- /dev/null +++ b/native/include/c_micro_sim.h @@ -0,0 +1,62 @@ +#pragma once +#include +#include +#include +#include +// #define MAX_ENERGY_TERMS 20 + +class Energy; // forward declaration; Energy is defined in c_energy.h +class MicroSim { +public: + MicroSim() {std::cout << "Instatiating MicroSim class; at " << this << "\n";}; + virtual ~MicroSim() {std::cout << "Killing MicroSim\n";}; + bool set_up; + + // From the mesh + int nx, ny, nz, n; + double dx, dy, dz; + double unit_length; + double *coordinates; + int *ngbs; + + // Arrays of material properties + double *spin; + double *Ms; + double *Ms_inv; + double *energy; + double *field; + int *pins; + + // Array with interactions + // void * interactions; + // int interactions_id[MAX_ENERGY_TERMS]; + // Should we use a vector of weak_ptr? : std::vector> interactions; + std::vector interactions; + + // Not necessary at least we re creating an array of void pointers: + // std::vector interactions_id; + + // Methods + void setup(int nx, int ny, int nz, double dx, double dy, double dz, + double unit_length, double *coordinates, int *ngbs, + double *spin, double *Ms, double *Ms_inv, + double *energy, double *field, int *pins + ); + + void add_interaction(Energy * interaction); + + // void print_interactions_id() { + // for(int i : interactions_id) std::cout << i << "\n"; + // for(auto i : interactions) std::cout << i << "\n"; + // } + + void compute_effective_field(double t); +}; + +enum EnergyTermIDs { + NONE = 0, + EXCHANGE, + DMI, + ZEEMAN, + UNIAXIAL_ANISOTROPY +}; diff --git a/fidimag/common/neb_method/nebm_geodesic_lib.h b/native/include/c_nebm_geodesic_lib.h similarity index 50% rename from fidimag/common/neb_method/nebm_geodesic_lib.h rename to native/include/c_nebm_geodesic_lib.h index 0de9c151..b3b03296 100644 --- a/fidimag/common/neb_method/nebm_geodesic_lib.h +++ b/native/include/c_nebm_geodesic_lib.h @@ -1,14 +1,13 @@ #include "math.h" -#define WIDE_PI 3.1415926535897932384626433832795L -double compute_geodesic_GreatCircle(double *restrict A, double *restrict B, +double compute_geodesic_GreatCircle(double * A, double * B, int n_dofs_image, - int *restrict material, + int * material, int n_dofs_image_material ); -double compute_geodesic_Vincenty(double *restrict A, double *restrict B, +double compute_geodesic_Vincenty(double * A, double * B, int n_dofs_image, - int *restrict material, + int * material, int n_dofs_image_material ); diff --git a/fidimag/common/neb_method/nebm_integrators.h b/native/include/c_nebm_integrators.h similarity index 51% rename from fidimag/common/neb_method/nebm_integrators.h rename to native/include/c_nebm_integrators.h index fac7e0f0..62fc34df 100644 --- a/fidimag/common/neb_method/nebm_integrators.h +++ b/native/include/c_nebm_integrators.h @@ -1,8 +1,8 @@ -double step_Verlet_C(double *restrict forces, - double *restrict forces_prev, - double *restrict velocities, - double *restrict velocities_new, - double *restrict y, +double step_Verlet_C(double * forces, + double * forces_prev, + double * velocities, + double * velocities_new, + double * y, double t, double h, double mass, diff --git a/native/include/c_nebm_lib.h b/native/include/c_nebm_lib.h new file mode 100644 index 00000000..612673cb --- /dev/null +++ b/native/include/c_nebm_lib.h @@ -0,0 +1,62 @@ +#include "math.h" + +void compute_tangents_C(double * ys, double * energy, + double * tangents, int image_num, int nodes + ); + +void compute_spring_force_C( + double * spring_force, + double * y, + double * tangents, + double * k, + int n_images, + int n_dofs_image, + double * distances + ); + +void compute_effective_force_C(double * G, + double * tangents, + double * gradientE, + double * spring_force, + int * climbing_image, + int n_images, + int n_dofs_image + ); + +void compute_image_distances(double * distances, + double * path_distances, + double * y, + int n_images, + int n_dofs_image, + double (* compute_distance)(double *, + double *, + int, + int *, + int), + int * material, + int n_dofs_image_material + ); + +void normalise_images_C(double * y, int n_images, int n_dofs_image); + +void normalise_spins_C(double * y, int n_images, int n_dofs_image); + +void project_images_C(double * vector, double * y, + int n_images, int n_dofs_image + ); + +void project_vector_C(double * vector, double * y, + int n_dofs_image + ); + +double compute_distance_cartesian(double * A, double * B, int n_dofs_image, + int * material, int n_dofs_image_material + ); + +void compute_dYdt_C(double * y, double * G, double * dYdt, + int * pins, + int n_images, int n_dofs_image); + +void compute_dYdt_nc_C(double * y, double * G, double * dYdt, + int * pins, + int n_images, int n_dofs_image); diff --git a/native/include/c_nebm_spherical_lib.h b/native/include/c_nebm_spherical_lib.h new file mode 100644 index 00000000..c18663c3 --- /dev/null +++ b/native/include/c_nebm_spherical_lib.h @@ -0,0 +1,10 @@ +#include "math.h" + +void normalise_spherical(double * a, int n); + +void normalise_images_spherical_C(double * y, int n_images, + int n_dofs_image); + +double compute_distance_spherical(double * A, double * B, int n, + int * material, int n_dofs_image_material + ); diff --git a/fidimag/common/dipolar/tensor_2dpbc.h b/native/include/c_tensor_2dpbc.h similarity index 97% rename from fidimag/common/dipolar/tensor_2dpbc.h rename to native/include/c_tensor_2dpbc.h index e2d41ae5..ca99f2cd 100644 --- a/fidimag/common/dipolar/tensor_2dpbc.h +++ b/native/include/c_tensor_2dpbc.h @@ -7,12 +7,12 @@ inline double RR(double x,double y,double z) inline double R(double x,double y,double z) {return sqrt(x*x+y*y+z*z);} -//------------------------------------------------------------------------------ +//------------------------------------------------------------------------------ //------------------------------------------------------------------------------ //about Nxxinf -inline double +inline double fxx(double x,double y,double z) { double t1,t2,t3,t5; @@ -49,7 +49,7 @@ inline double fzz(double x,double y,double z){ } -inline double +inline double Nzzinf(double x,double y,double z,double X0,double Y0) { return fzz(x+X0,y+Y0,z)+fzz(x-X0,y-Y0,z)-fzz(x+X0,y-Y0,z)-fzz(x-X0,y+Y0,z); @@ -97,7 +97,7 @@ inline double Nyzinf(double x,double y,double z,double X0,double Y0){ //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ //about Nxxdipole -inline double +inline double DemagNxxDipolar(double x,double y,double z) { double t1 = x*x; @@ -109,7 +109,7 @@ double t8 = sqrt(t6); return -(2.0 * t1 - t3 - t4) / (t8 * t7); } -inline double +inline double DemagNxyDipolar(double x,double y,double z) { double t6 = RR(x,y,z); @@ -125,4 +125,3 @@ double DemagTensorAsymptotic(enum Type_Nij comp,double x,double y,double z,doubl double DemagTensorDipolar(enum Type_Nij comp,double x,double y,double z); double DemagTensorInfinite(enum Type_Nij comp,double x,double y,double z,double X0,double Y0); - diff --git a/native/include/c_vectormath.h b/native/include/c_vectormath.h new file mode 100644 index 00000000..5a28a600 --- /dev/null +++ b/native/include/c_vectormath.h @@ -0,0 +1,15 @@ +#pragma once +#include + + +void cross(double *output, double * A, double * B); +double dot(double *A, double *B, int n); + +double compute_norm(double * a, int n); + +void normalise(double * a, int n); +void normalise(double *m, int *pins, int n); + +double cross_x(double a0, double a1, double a2, double b0, double b1, double b2); +double cross_y(double a0, double a1, double a2, double b0, double b1, double b2); +double cross_z(double a0, double a1, double a2, double b0, double b1, double b2); diff --git a/fidimag/micro/lib/baryakhtar/baryakhtar_clib.h b/native/include/m_baryakhtar_clib.h similarity index 65% rename from fidimag/micro/lib/baryakhtar/baryakhtar_clib.h rename to native/include/m_baryakhtar_clib.h index 99c7fec3..cec829ea 100644 --- a/fidimag/micro/lib/baryakhtar/baryakhtar_clib.h +++ b/native/include/m_baryakhtar_clib.h @@ -1,14 +1,7 @@ #include -//#include - -#define WIDE_PI 3.1415926535897932384626433832795L #define MU0 1.25663706143591728850e-6 -inline double cross_x(double a0, double a1, double a2, double b0, double b1, double b2) { return a1*b2 - a2*b1; } -inline double cross_y(double a0, double a1, double a2, double b0, double b1, double b2) { return a2*b0 - a0*b2; } -inline double cross_z(double a0, double a1, double a2, double b0, double b1, double b2) { return a0*b1 - a1*b0; } - void compute_laplace_m(double *m, double *field, double *Ms, double dx, double dy, double dz, int nx, int ny, int nz); diff --git a/native/include/m_clib.h b/native/include/m_clib.h new file mode 100644 index 00000000..49e00aa7 --- /dev/null +++ b/native/include/m_clib.h @@ -0,0 +1,22 @@ +#include + +#define MU0 1.25663706143591728850e-6 +#define MU0_INV 795774.71545947669074 + +void compute_exch_field_micro(double * m, double * field, double * energy, double * Ms_inv, + double A, double dx, double dy, double dz, int n, int *ngbs); + +void dmi_field(double * m, double * field, + double * energy, double * Ms_inv, + double * D, int n_DMIs, + double *dmi_vector, + double dx, double dy, double dz, int n, int *ngbs); + +void compute_exch_field_rkky_micro(double *m, double *field, double *energy, double *Ms_inv, + double sigma, int nx, double ny, double nz, int z_bottom, int z_top); + +void compute_uniaxial_anis(double * m, double * field, double * energy, double * Ms_inv, + double * Ku, double * axis, int nx, int ny, int nz); + +double skyrmion_number(double * spin, double * charge, + int nx, int ny, int nz, int * ngbs); diff --git a/native/include/m_driver.h b/native/include/m_driver.h new file mode 100644 index 00000000..ccdea839 --- /dev/null +++ b/native/include/m_driver.h @@ -0,0 +1,126 @@ +#pragma once +#include +#include +#include "c_micro_sim.h" + +// TODO: We should move this to a separate file containing all the structs +struct TestParams { + double foo; +}; + +enum IntegratorID { + RK4 +}; + +// We will set the extra parameters for the integrator as a struct +// The integrator allows m, dmdt, time, dt and extra parameters +struct LLG_params { + double gamma; + double * alpha; + double * field; + int * pins; +}; + +// In case we want to use conditional statements for LLG eq: +// enum LLGTerms { +// PRECESSION, +// DAMPING, +// STT_ZHANGLI, +// STT_SLONCZEWSKI, +// STT_FIELDLIKE +// }; + +template +class Integrator; // forward declaration + +// Base class for the drivers -> just declare LLG driver for now +class MicroLLGDriver { +public: + MicroLLGDriver() {std::cout << "Instatiating LLG Driver" << std::endl;}; + ~MicroLLGDriver() {std::cout << "Killing LLG Driver\n";}; + + // Testing constructor + // MicroLLGDriver(MicroSim * sim, double * alpha, double gamma, + // double t, double dt) { + // sim = sim; + // alpha = alpha; + // gamma = gamma; + // t = t; + // dt = dt; + // } + + // TODO: declare all remaining values from the full form LLG eq + // Parameters are set in the LLG_params struct + LLG_params * llg_params; + double t; + double dt; + MicroSim * sim; + // This creates a pointer to the base class Integrator; with add_integrator + // the virtual functions are override by corresp. derived class functions + // of the chosen integrator + Integrator * integrator; + std::vector dmdt; // N len vector, we could also use an array + + // Will get the parameters from a simulation class + void setup(MicroSim * sim, + double * alpha, double gamma, + // double * zl_jx, double * zl_jy, double * zl_jz, double zl_p, double zl_beta, double zl_u0, + // double * cpp_p, double * cpp_aJ, double * cpp_beta, + double t, double dt); + void add_integrator(IntegratorID integrator_id); + void run_until(double t_final); + void single_integrator_step(); + + // void compute_RHS(double * m, std::vector& dmdt, + // unsigned int n, double t); +}; + +// ---------------------------------------------------------------------------- +// TODO: Integrators should be moved to a separate file +// +// Abstract class for the integrator +// The template refers to the type of the pointer to pass extra parameters to +// the integrator, e.g. the LLG Driver uses a pointer to the LLG_params struct +template +class Integrator { +public: + // Integrator() = {}; + // We should probably better use a constructor: + // Accepts N variables (e.g. 3 * n_spins) and a pointer to a function + virtual void setup(unsigned int N, double dt) = 0; + + // Templated pointer for parameters, e.g. for LLG we need eff field, alpha, etc + virtual void integration_step(void (*f) (double * m, std::vector& dmdt, + unsigned int n, double t, T * params), + double * m, std::vector& dmdt, + T * params) = 0; + + IntegratorID integrator_id; + + unsigned int step_n; + unsigned int N; + double t; + double dt; + + std::vector integratorData; +}; + +// Integrators should be independent of any driver/sim class +template +class IntegratorRK4: public Integrator { +public: + IntegratorRK4(size_t N) { + // Know we need 4xN for RK4 + this->integratorData.resize(4 * N); + std::cout << "Instatiating RK4" << std::endl; + }; + + virtual ~IntegratorRK4() {std::cout << "Killing RK4 integrator\n";}; + + // std::vector rksteps; // N len vector, we could also use an array + + void setup(unsigned int N, double dt); + void integration_step(void (*f) (double * m, std::vector& dmdt, + unsigned int N, double t, T * params), + double * m, std::vector& dmdt, T * params); +}; diff --git a/native/include/m_driver_test.h b/native/include/m_driver_test.h new file mode 100644 index 00000000..237f534c --- /dev/null +++ b/native/include/m_driver_test.h @@ -0,0 +1,4 @@ +#pragma once +#include "m_driver.h" + +void test_integrator_RK4 (void); diff --git a/fidimag/atomistic/lib/anis.c b/native/src/a_anis.cpp similarity index 79% rename from fidimag/atomistic/lib/anis.c rename to native/src/a_anis.cpp index 73217e7f..be12ca0a 100644 --- a/fidimag/atomistic/lib/anis.c +++ b/native/src/a_anis.cpp @@ -1,10 +1,10 @@ -#include "clib.h" +#include "a_clib.h" -void compute_anis(double *restrict spin, double *restrict field, - double *restrict mu_s_inv, - double *restrict energy, - double *restrict Ku, double *restrict axis, int n) { +void compute_anis(double * spin, double * field, + double * mu_s_inv, + double * energy, + double * Ku, double * axis, int n) { /* Remember that the magnetisation order is * mx1, my1, mz1, mx2, my2, mz2, mx3,... @@ -40,10 +40,10 @@ void compute_anis(double *restrict spin, double *restrict field, } -void compute_anis_cubic(double *restrict spin, double *restrict field, - double *restrict mu_s_inv, - double *restrict energy, - double *restrict Kc, int n) { +void compute_anis_cubic(double * spin, double * field, + double * mu_s_inv, + double * energy, + double * Kc, int n) { /* Remember that the magnetisation order is * mx1, my1, mz1, mx2, my2, mz2, mx3,... diff --git a/fidimag/atomistic/lib/demag_full.c b/native/src/a_demag_full.cpp similarity index 94% rename from fidimag/atomistic/lib/demag_full.c rename to native/src/a_demag_full.cpp index 69432e78..f387588b 100644 --- a/fidimag/atomistic/lib/demag_full.c +++ b/native/src/a_demag_full.cpp @@ -1,9 +1,9 @@ -#include "clib.h" +#include "a_clib.h" #include "math.h" #include "stdlib.h" -void demag_full(double *restrict spin, double *restrict field, double *restrict energy, double *restrict coords, - double *restrict mu_s, double *restrict mu_s_scale, int n) { +void demag_full(double * spin, double * field, double * energy, double * coords, + double * mu_s, double * mu_s_scale, int n) { /* Full calculation of the Demagnetising field for atomistic systems * The main idea is to iterate through every lattice site, and sum diff --git a/fidimag/atomistic/lib/dmi.c b/native/src/a_dmi.cpp similarity index 92% rename from fidimag/atomistic/lib/dmi.c rename to native/src/a_dmi.cpp index b53b8a2a..d98892f6 100644 --- a/fidimag/atomistic/lib/dmi.c +++ b/native/src/a_dmi.cpp @@ -1,6 +1,7 @@ -#include "clib.h" +#include "a_clib.h" #include "math.h" #include "stdlib.h" +#include "c_vectormath.h" /* NOTES about the neighbours array in the arguments: @@ -27,10 +28,10 @@ boundaries */ -void dmi_field_bulk(double *restrict spin, double *restrict field, - double *restrict mu_s_inv, - double *restrict energy, double *restrict _D, - int *restrict ngbs, int nxyz, int n_ngbs) { +void dmi_field_bulk(double * spin, double * field, + double * mu_s_inv, + double * energy, double * _D, + int * ngbs, int nxyz, int n_ngbs) { /* Bulk DMI field and energy computation * @@ -119,12 +120,12 @@ void dmi_field_bulk(double *restrict spin, double *restrict field, } } -void dmi_field_interfacial_atomistic(double *restrict spin, double *restrict field, - double *restrict mu_s_inv, - double *restrict energy, - double D, int *restrict ngbs, int n, +void dmi_field_interfacial_atomistic(double * spin, double * field, + double * mu_s_inv, + double * energy, + double D, int * ngbs, int n, int n_ngbs, int n_ngbs_dmi, - double *restrict DMI_vec) { + double * DMI_vec) { /* Interfacial DMI field and energy computation * @@ -233,7 +234,7 @@ inline double single_energy_z(double D, double Si[3], double Sj[3]){ return tz; } -double dmi_energy(double *restrict spin, double D, int nx, int ny, int nz, int xperiodic, int yperiodic) { +double dmi_energy(double * spin, double D, int nx, int ny, int nz, int xperiodic, int yperiodic) { int nyz = ny * nz; int n1 = nx * nyz, n2 = 2 * n1; diff --git a/fidimag/atomistic/lib/exch.c b/native/src/a_exch.cpp similarity index 88% rename from fidimag/atomistic/lib/exch.c rename to native/src/a_exch.cpp index 81d65210..f38386ce 100644 --- a/fidimag/atomistic/lib/exch.c +++ b/native/src/a_exch.cpp @@ -1,4 +1,4 @@ -#include "clib.h" +#include "a_clib.h" /* compute the effective exchange field at site i @@ -35,11 +35,11 @@ The ngbs array also gives the correct indexes for the spins at periodic boundaries */ -void compute_exch_field(double *restrict spin, double *restrict field, - double *restrict mu_s_inv, - double *restrict energy, +void compute_exch_field(double * spin, double * field, + double * mu_s_inv, + double * energy, double Jx, double Jy, double Jz, - int *restrict ngbs, int n, int n_ngbs) { + int * ngbs, int n, int n_ngbs) { #pragma omp parallel for for (int i = 0; i < n; i++) { @@ -76,7 +76,7 @@ void compute_exch_field(double *restrict spin, double *restrict field, -double compute_exch_energy(double *restrict spin, +double compute_exch_energy(double * spin, double Jx, double Jy, double Jz, int nx, int ny, int nz, int xperiodic, int yperiodic) { @@ -145,10 +145,10 @@ double compute_exch_energy(double *restrict spin, Note that the pair only run once for each pair. */ -void compute_exch_field_spatial(double *restrict spin, double *restrict field, - double *restrict mu_s_inv, - double *restrict energy, - double *restrict J, int *restrict ngbs, +void compute_exch_field_spatial(double * spin, double * field, + double * mu_s_inv, + double * energy, + double * J, int * ngbs, int n, int n_ngbs) { #pragma omp parallel for @@ -203,12 +203,12 @@ void compute_exch_field_spatial(double *restrict spin, double *restrict field, * Thus, we can locate ngbs from cols 0-5, 6-11, 12, 23, ... etc * */ -void compute_full_exch_field(double *restrict spin, double *restrict field, - double *restrict mu_s_inv, - double *restrict energy, +void compute_full_exch_field(double * spin, double * field, + double * mu_s_inv, + double * energy, double J[9], int *ngbs, int n, int n_ngbs, - int n_shells, int *restrict n_ngbs_shell, - int *restrict sum_ngbs_shell + int n_shells, int * n_ngbs_shell, + int * sum_ngbs_shell ) { #pragma omp parallel for diff --git a/fidimag/atomistic/lib/mc.c b/native/src/a_mc.cpp similarity index 97% rename from fidimag/atomistic/lib/mc.c rename to native/src/a_mc.cpp index e19a37a9..a7395906 100644 --- a/fidimag/atomistic/lib/mc.c +++ b/native/src/a_mc.cpp @@ -1,5 +1,5 @@ -#include "clib.h" -#include "fidimag_random.h" +#include "a_clib.h" +#include "a_random.h" inline double dot(double a[3], double b[3]) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; @@ -75,7 +75,7 @@ inline double cubic_energy_site(double *m, double Kc) { return -Kc * (mx2 * mx2 + my2 * my2 + mz2 * mz2); } -double compute_deltaE_anisotropy(double *restrict spin, double *restrict new_spin, double *restrict h, +double compute_deltaE_anisotropy(double * spin, double * new_spin, double * h, double Kc, int i, int new_i) { double energy1 = -dot(&spin[3 * i], &h[3 * i]); // zeeman energy diff --git a/fidimag/atomistic/lib/fidimag_random.c b/native/src/a_random.cpp similarity index 97% rename from fidimag/atomistic/lib/fidimag_random.c rename to native/src/a_random.cpp index dda39853..3bad4832 100644 --- a/fidimag/atomistic/lib/fidimag_random.c +++ b/native/src/a_random.cpp @@ -1,16 +1,17 @@ #include // rand(), srand() #include #include -#include "fidimag_random.h" +#include "a_random.h" +#include "c_constants.h" mt19937_state *create_mt19937_state(void) { - + mt19937_state *state = (mt19937_state*) malloc(sizeof(mt19937_state)); state->seed = 1; state->index_t = 0; state->matrix[0] = 0; state->matrix[1] = 0x9908b0dfU; - + return state; } @@ -20,13 +21,13 @@ void finalize_mt19937_state(mt19937_state *state) { void initial_rng_mt19973(mt19937_state *state, int seed) { - + if (seed<0){ state->seed = (unsigned int) time(NULL); }else{ state->seed = seed; } - + state->MT[0] = state->seed & 0xFFFFFFFFU; for (int i = 1; i < 624; i++) { state->MT[i] = (state->MT[i - 1] ^ (state->MT[i - 1] >> 30)) + i; @@ -37,21 +38,21 @@ void initial_rng_mt19973(mt19937_state *state, int seed) { //return a integer in [0, MT19973_RAND_MAX] inline unsigned int rand_int(mt19937_state *state) { - + unsigned int x; int index_t = state->index_t; - + x = (state->MT[index_t] & 0x1U) + (state->MT[(index_t + 1) % 624] & 0xFFFFFFFEU); state->MT[index_t] = (state->MT[(index_t + 397) % 624] ^ (x >> 1))^ state->matrix[x & 1]; - + x = state->MT[index_t]; x ^= (x >> 11); x ^= (x << 7) & 0x9d2c5680U; x ^= (x << 15) & 0xefc60000U; x ^= (x >> 18); - + state->index_t = (index_t + 1) % 624; - + return x; } @@ -135,7 +136,7 @@ void random_integer_vector(mt19937_state *state, int *ids, int n) { void uniform_random_sphere(mt19937_state *state, double *spin, int n){ for (int i=0;i0) continue; - - - mm = 1.0 / sqrt(m[i] * m[i] + m[j] * m[j] + m[k] * m[k]); - m[i] *= mm; - m[j] *= mm; - m[k] *= mm; - - } -} diff --git a/fidimag/atomistic/lib/util.c b/native/src/a_util.cpp similarity index 91% rename from fidimag/atomistic/lib/util.c rename to native/src/a_util.cpp index d38d1da5..eba806f7 100644 --- a/fidimag/atomistic/lib/util.c +++ b/native/src/a_util.cpp @@ -1,6 +1,10 @@ -#include "clib.h" -#include "complex.h" -#include "math.h" +#include "a_clib.h" +#include "c_vectormath.h" +#include +#include +#include "c_constants.h" + +using namespace std::complex_literals; // compute the S \cdot (S_i \times S_j) inline double volume(double S[3], double Si[3], double Sj[3]) { @@ -10,8 +14,8 @@ inline double volume(double S[3], double Si[3], double Sj[3]) { return tx + ty + tz; } -double skyrmion_number(double *restrict spin, double *restrict charge, int nx, int ny, int nz, - int *restrict ngbs, int n_ngbs) { +double skyrmion_number(double * spin, double * charge, int nx, int ny, int nz, + int * ngbs, int n_ngbs) { /* Calculation of the "Skyrmion number" Q for a two dimensional discrete * spin lattice in the x-y plane (also known @@ -31,7 +35,7 @@ double skyrmion_number(double *restrict spin, double *restrict charge, int nx, i * lattice site. The array is like: * * __ indexes beyond nearest neighbours - * | + * | * [ 0-x, 0+x, 0-y, 0+y, 0-z, 0+z, ... 1-x, 1+x, 1-y, ... ] * i=0 i=n_ngbs ... * @@ -151,7 +155,7 @@ double skyrmion_number(double *restrict spin, double *restrict charge, int nx, i charge[i] += volume(S, S_i, S_j); /* Scale the chirality quantity */ - charge[i] /= (8 * WIDE_PI); + charge[i] /= (8 * M_PIl); /* We use the sum to output the total spin chirality * or skyrmion number */ @@ -161,20 +165,6 @@ double skyrmion_number(double *restrict spin, double *restrict charge, int nx, i return sum; } -double dot(double *a, double *b) { - /* Dot product for vectors of 3 components, given by arrays a and b. If - * pointers are passed, it uses the first three components starting from - * that place in memory, i.e. if we pass &a[2], where a = {0, 1, 2, 3, 4}, - * the dot product uses {2, 3, 4} as a vector (same for b) - */ - double dp = 0; - int i; - - for (i = 0; i < 3; i++) - dp += a[i] * b[i]; - return dp; -} - double compute_BergLuscher_angle(double *s1, double *s2, double *s3) { /* Compute the spherical angle given by the neighbouring spin 3-vectors s1, @@ -182,11 +172,9 @@ double compute_BergLuscher_angle(double *s1, double *s2, double *s3) { * -dot- function). The spherical angle Omega, defined by the * vectors in a unit sphere, is computed using the imaginary exponential * defined by Berg and Luscher [Nucl Phys B 190, 412 (1981)]: - * exp (i sigma Omega) = 1 + s1 * s2 + s2 * s3 + s3 * s1 + i s1 * (s2 X s3) * ------------------------------------------------ * 2 (1 + s1 * s2) (1 + s2 * s3) (1 + s3 * s1) - * The denominator is a normalisation factor which we call rho, and i * stands for an imaginary number in the numerator. The factor sigma is an * orientation given by sign(s1 * s2 X s3), however, we do not use it since @@ -194,9 +182,7 @@ double compute_BergLuscher_angle(double *s1, double *s2, double *s3) { * spins, as pointed out by Yin et al. [PRB 93, 174403 (2016)]. * * Therefore we use a complex logarithm: - * clog( r * exp(i theta) ) = r + i theta - * to calculate the angle Omega, since the clog is well defined in the * [-PI, PI] range, giving the correct sign for the topological number (we * could also use the arcsin when decomposing the exp). @@ -206,21 +192,19 @@ double compute_BergLuscher_angle(double *s1, double *s2, double *s3) { */ double rho; - double complex exp; + std::complex exp; double crossp[3]; - crossp[0] = s2[1] * s3[2] - s2[2] * s3[1]; - crossp[1] = s2[2] * s3[0] - s2[0] * s3[2]; - crossp[2] = s2[0] * s3[1] - s2[1] * s3[0]; + cross(crossp, s2, s3); - rho = sqrt(2 * (1 + dot(&s1[0], &s2[0])) * (1 + dot(&s2[0], &s3[0])) * - (1 + dot(&s3[0], &s1[0]))); + rho = sqrt(2 * (1 + dot(&s1[0], &s2[0], 3)) * (1 + dot(&s2[0], &s3[0], 3)) * + (1 + dot(&s3[0], &s1[0], 3))); - exp = (1 + dot(&s1[0], &s2[0]) + dot(&s2[0], &s3[0]) + dot(&s3[0], &s1[0]) + - I * dot(&s1[0], &crossp[0])) / + exp = (1 + dot(&s1[0], &s2[0], 3) + dot(&s2[0], &s3[0], 3) + dot(&s3[0], &s1[0], 3) + + 1i * dot(&s1[0], &crossp[0], 3)) / rho; - return 2 * cimagl(clog(exp)) / (4 * WIDE_PI); + return 2 * std::imag(std::log(exp)) / (4 * M_PIl); } double skyrmion_number_BergLuscher(double *spin, double *charge, int nx, int ny, @@ -240,8 +224,8 @@ double skyrmion_number_BergLuscher(double *spin, double *charge, int nx, int ny, * site. The first 6 elements are the NEAREST neighbours. For cuboid meshes, * the array is like: * - * __ indexes beyond nearest ngbs - * | + * __ indexes beyond nearest ngbs + * | * NN: j =0 j=1 ... | j=0 j=1 ... * [ 0-x, 0+x, 0-y, 0+y, 0-z, 0+z, ... 1-x, 1+x, 1-y, ... ] * spin: i=n_ngbs * 0 i=n_ngbs * 1 diff --git a/fidimag/common/dipolar/demag.c b/native/src/c_demag.cpp similarity index 67% rename from fidimag/common/dipolar/demag.c rename to native/src/c_demag.cpp index 1c908ae2..67202f4e 100644 --- a/fidimag/common/dipolar/demag.c +++ b/native/src/c_demag.cpp @@ -1,8 +1,8 @@ #include #include #include -#include "dipolar.h" -#include "demagcoef.h" +#include "c_dipolar.h" +#include "c_demagcoef.h" inline double Nxxdipole(double x, double y, double z) { @@ -165,7 +165,7 @@ void compute_demag_tensors(fft_demag_plan *plan) { //used for debug -void print_r(char *str, double *restrict x, int n) { +void print_r(char *str, double * x, int n) { int i; printf("%s:\n", str); for (i = 0; i < n; i++) { @@ -175,7 +175,7 @@ void print_r(char *str, double *restrict x, int n) { } -void print_c(char *str, fftw_complex *restrict x, int n) { +void print_c(char *str, fftw_complex * x, int n) { int i; printf("%s\n", str); for (i = 0; i < n; i++) { @@ -192,7 +192,7 @@ fft_demag_plan *create_plan(void) { return plan; } -void init_plan(fft_demag_plan *restrict plan, double dx, double dy, +void init_plan(fft_demag_plan * plan, double dx, double dy, double dz, int nx, int ny, int nz) { //plan->mu_s = mu_s; @@ -217,7 +217,7 @@ void init_plan(fft_demag_plan *restrict plan, double dx, double dy, plan->total_length = plan->lenx * plan->leny * plan->lenz; int size1 = plan->total_length * sizeof(double); - int size2 = plan->total_length * sizeof(fftw_complex); + int size2 = plan->total_length * sizeof(std::complex); plan->tensor_xx = (double *) fftw_malloc(size1); plan->tensor_yy = (double *) fftw_malloc(size1); @@ -234,36 +234,36 @@ void init_plan(fft_demag_plan *restrict plan, double dx, double dy, plan->hy = (double *) fftw_malloc(size1); plan->hz = (double *) fftw_malloc(size1); - plan->Nxx = (fftw_complex *) fftw_malloc(size2); - plan->Nyy = (fftw_complex *) fftw_malloc(size2); - plan->Nzz = (fftw_complex *) fftw_malloc(size2); - plan->Nxy = (fftw_complex *) fftw_malloc(size2); - plan->Nxz = (fftw_complex *) fftw_malloc(size2); - plan->Nyz = (fftw_complex *) fftw_malloc(size2); + plan->Nxx = (std::complex *) fftw_malloc(size2); + plan->Nyy = (std::complex *) fftw_malloc(size2); + plan->Nzz = (std::complex *) fftw_malloc(size2); + plan->Nxy = (std::complex *) fftw_malloc(size2); + plan->Nxz = (std::complex *) fftw_malloc(size2); + plan->Nyz = (std::complex *) fftw_malloc(size2); - plan->Mx = (fftw_complex *) fftw_malloc(size2); - plan->My = (fftw_complex *) fftw_malloc(size2); - plan->Mz = (fftw_complex *) fftw_malloc(size2); - plan->Hx = (fftw_complex *) fftw_malloc(size2); - plan->Hy = (fftw_complex *) fftw_malloc(size2); - plan->Hz = (fftw_complex *) fftw_malloc(size2); + plan->Mx = (std::complex *) fftw_malloc(size2); + plan->My = (std::complex *) fftw_malloc(size2); + plan->Mz = (std::complex *) fftw_malloc(size2); + plan->Hx = (std::complex *) fftw_malloc(size2); + plan->Hy = (std::complex *) fftw_malloc(size2); + plan->Hz = (std::complex *) fftw_malloc(size2); } -void create_fftw_plan(fft_demag_plan *restrict plan) { +void create_fftw_plan(fft_demag_plan * plan) { plan->tensor_plan = fftw_plan_dft_r2c_3d(plan->lenz, plan->leny, - plan->lenx, plan->tensor_xx, plan->Nxx, - FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); - + plan->lenx, plan->tensor_xx, reinterpret_cast(plan->Nxx), + FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); + plan->m_plan = fftw_plan_dft_r2c_3d(plan->lenz, plan->leny, plan->lenx, - plan->mx, plan->Mx, FFTW_MEASURE); - + plan->mx, reinterpret_cast(plan->Mx), FFTW_MEASURE); + plan->h_plan = fftw_plan_dft_c2r_3d(plan->lenz, plan->leny, plan->lenx, - plan->Hx, plan->hx, FFTW_MEASURE | FFTW_DESTROY_INPUT); - + reinterpret_cast(plan->Hx), plan->hx, FFTW_MEASURE | FFTW_DESTROY_INPUT); + for (int i = 0; i < plan->total_length; i++) { plan->Nxx[i] = 0; plan->Nyy[i] = 0; @@ -282,12 +282,12 @@ void create_fftw_plan(fft_demag_plan *restrict plan) { } - fftw_execute_dft_r2c(plan->tensor_plan, plan->tensor_xx, plan->Nxx); - fftw_execute_dft_r2c(plan->tensor_plan, plan->tensor_yy, plan->Nyy); - fftw_execute_dft_r2c(plan->tensor_plan, plan->tensor_zz, plan->Nzz); - fftw_execute_dft_r2c(plan->tensor_plan, plan->tensor_xy, plan->Nxy); - fftw_execute_dft_r2c(plan->tensor_plan, plan->tensor_xz, plan->Nxz); - fftw_execute_dft_r2c(plan->tensor_plan, plan->tensor_yz, plan->Nyz); + fftw_execute_dft_r2c(plan->tensor_plan, plan->tensor_xx, reinterpret_cast(plan->Nxx)); + fftw_execute_dft_r2c(plan->tensor_plan, plan->tensor_yy, reinterpret_cast(plan->Nyy)); + fftw_execute_dft_r2c(plan->tensor_plan, plan->tensor_zz, reinterpret_cast(plan->Nzz)); + fftw_execute_dft_r2c(plan->tensor_plan, plan->tensor_xy, reinterpret_cast(plan->Nxy)); + fftw_execute_dft_r2c(plan->tensor_plan, plan->tensor_xz, reinterpret_cast(plan->Nxz)); + fftw_execute_dft_r2c(plan->tensor_plan, plan->tensor_yz, reinterpret_cast(plan->Nyz)); fftw_destroy_plan(plan->tensor_plan); } @@ -298,7 +298,7 @@ void create_fftw_plan(fft_demag_plan *restrict plan) { //The computed results doesn't consider the coefficient of \frac{\mu_0}{4 \pi}, the //reason is in future we can use the following code directly for continuum case -void compute_fields(fft_demag_plan *restrict plan, double *restrict spin, double *restrict mu_s, double *restrict field) { +void compute_fields(fft_demag_plan * plan, double * spin, double * mu_s, double * field) { int i, j, k, id1, id2; @@ -333,25 +333,24 @@ void compute_fields(fft_demag_plan *restrict plan, double *restrict spin, double //print_r("plan->mx", plan->mx, plan->total_length); - fftw_execute_dft_r2c(plan->m_plan, plan->mx, plan->Mx); - fftw_execute_dft_r2c(plan->m_plan, plan->my, plan->My); - fftw_execute_dft_r2c(plan->m_plan, plan->mz, plan->Mz); + fftw_execute_dft_r2c(plan->m_plan, plan->mx, reinterpret_cast(plan->Mx)); + fftw_execute_dft_r2c(plan->m_plan, plan->my, reinterpret_cast(plan->My)); + fftw_execute_dft_r2c(plan->m_plan, plan->mz, reinterpret_cast(plan->Mz)); //print_c("plan->Mx", plan->Mx, plan->total_length); - - fftw_complex *Nxx = plan->Nxx; - fftw_complex *Nyy = plan->Nyy; - fftw_complex *Nzz = plan->Nzz; - fftw_complex *Nxy = plan->Nxy; - fftw_complex *Nxz = plan->Nxz; - fftw_complex *Nyz = plan->Nyz; - - fftw_complex *Mx = plan->Mx; - fftw_complex *My = plan->My; - fftw_complex *Mz = plan->Mz; - fftw_complex *Hx = plan->Hx; - fftw_complex *Hy = plan->Hy; - fftw_complex *Hz = plan->Hz; + std::complex *Nxx = plan->Nxx; + std::complex *Nyy = plan->Nyy; + std::complex *Nzz = plan->Nzz; + std::complex *Nxy = plan->Nxy; + std::complex *Nxz = plan->Nxz; + std::complex *Nyz = plan->Nyz; + + std::complex *Mx = plan->Mx; + std::complex *My = plan->My; + std::complex *Mz = plan->Mz; + std::complex *Hx = plan->Hx; + std::complex *Hy = plan->Hy; + std::complex *Hz = plan->Hz; //print_c("Mx", Mx, plan->total_length); @@ -363,9 +362,9 @@ void compute_fields(fft_demag_plan *restrict plan, double *restrict spin, double //print_c("Hx", Hx, plan->total_length); - fftw_execute_dft_c2r(plan->h_plan, plan->Hx, plan->hx); - fftw_execute_dft_c2r(plan->h_plan, plan->Hy, plan->hy); - fftw_execute_dft_c2r(plan->h_plan, plan->Hz, plan->hz); + fftw_execute_dft_c2r(plan->h_plan, reinterpret_cast(plan->Hx), plan->hx); + fftw_execute_dft_c2r(plan->h_plan, reinterpret_cast(plan->Hy), plan->hy); + fftw_execute_dft_c2r(plan->h_plan, reinterpret_cast(plan->Hz), plan->hz); //print_r("hx", plan->hx, plan->total_length); //print_r("hy", plan->hy, plan->total_length); //print_r("hz", plan->hz, plan->total_length); @@ -388,78 +387,78 @@ void compute_fields(fft_demag_plan *restrict plan, double *restrict spin, double } //only used for debug -void exact_compute(fft_demag_plan *restrict plan, double *restrict spin, double *restrict mu_s, double *restrict field) { - int i, j, k, index; - int ip, jp, kp, idf, ids; - int nx = plan->nx; - int ny = plan->ny; - int nz = plan->nz; - int nxy = nx * ny; - - //int lenx = plan->lenx; - int lenx = plan->lenx; - int leny = plan->leny; - int lenz = plan->lenz; - int lenxy = lenx * leny; - - double *Nxx = plan->tensor_xx; - double *Nyy = plan->tensor_yy; - double *Nzz = plan->tensor_zz; - double *Nxy = plan->tensor_xy; - double *Nxz = plan->tensor_xz; - double *Nyz = plan->tensor_yz; - +void exact_compute(fft_demag_plan * plan, double * spin, double * mu_s, double * field) { + int i, j, k, index; + int ip, jp, kp, idf, ids; + int nx = plan->nx; + int ny = plan->ny; + int nz = plan->nz; + int nxy = nx * ny; + + //int lenx = plan->lenx; + int lenx = plan->lenx; + int leny = plan->leny; + int lenz = plan->lenz; + int lenxy = lenx * leny; + + double *Nxx = plan->tensor_xx; + double *Nyy = plan->tensor_yy; + double *Nzz = plan->tensor_zz; + double *Nxy = plan->tensor_xy; + double *Nxz = plan->tensor_xz; + double *Nyz = plan->tensor_yz; + + + for (k = 0; k < nz; k++) { + for (j = 0; j < ny; j++) { + for (i = 0; i < nx; i++) { + idf = nxy * k + nx * j + i; - for (k = 0; k < nz; k++) { - for (j = 0; j < ny; j++) { - for (i = 0; i < nx; i++) { - idf = nxy * k + nx * j + i; - - field[3*idf] = 0; - field[3*idf+1] = 0; - field[3*idf+2] = 0; - - for (kp = 0; kp < nz; kp++) { - for (jp = 0; jp < ny; jp++) { - for (ip = 0; ip < nx; ip++) { - ids = nxy * kp + nx * jp + ip; - - index = ip-i>=0 ? ip - i: ip - i + lenx; - index += jp-j >=0 ? (jp - j)*lenx: (jp - j + leny)*lenx; - index += kp-k >=0 ? (kp - k)*lenxy: (kp - k + lenz)*lenxy; - - field[3*idf] += (Nxx[index] * spin[3*ids] + Nxy[index] - * spin[3*ids+1] + Nxz[index] * spin[3*ids+2])*mu_s[ids]; - field[3*idf+1] += (Nxy[index] * spin[3*ids] + Nyy[index] - * spin[3*ids+1] + Nyz[index] * spin[3*ids+2])*mu_s[ids]; - field[3*idf+2] += (Nxz[index] * spin[3*ids] + Nyz[index] - * spin[3*ids+1] + Nzz[index] * spin[3*ids+2])*mu_s[ids]; - } - } - } - - field[3*idf] *= (-1); - field[3*idf+1] *= (-1); - field[3*idf+2] *= (-1); - } - } + field[3*idf] = 0; + field[3*idf+1] = 0; + field[3*idf+2] = 0; + + for (kp = 0; kp < nz; kp++) { + for (jp = 0; jp < ny; jp++) { + for (ip = 0; ip < nx; ip++) { + ids = nxy * kp + nx * jp + ip; + + index = ip-i>=0 ? ip - i: ip - i + lenx; + index += jp-j >=0 ? (jp - j)*lenx: (jp - j + leny)*lenx; + index += kp-k >=0 ? (kp - k)*lenxy: (kp - k + lenz)*lenxy; + + field[3*idf] += (Nxx[index] * spin[3*ids] + Nxy[index] + * spin[3*ids+1] + Nxz[index] * spin[3*ids+2])*mu_s[ids]; + field[3*idf+1] += (Nxy[index] * spin[3*ids] + Nyy[index] + * spin[3*ids+1] + Nyz[index] * spin[3*ids+2])*mu_s[ids]; + field[3*idf+2] += (Nxz[index] * spin[3*ids] + Nyz[index] + * spin[3*ids+1] + Nzz[index] * spin[3*ids+2])*mu_s[ids]; + } + } } - + + field[3*idf] *= (-1); + field[3*idf+1] *= (-1); + field[3*idf+2] *= (-1); + } + } + } + } -double compute_demag_energy(fft_demag_plan *restrict plan, - double *restrict spin, - double *restrict mu_s, - double *restrict field, - double *restrict energy +double compute_demag_energy(fft_demag_plan * plan, + double * spin, + double * mu_s, + double * field, + double * energy ) { - + int i,j; - + int nxyz = plan->nx * plan->ny * plan->nz; - + double total_energy = 0; - + for (i = 0; i < nxyz; i++) { j = 3*i; energy[i] = -0.5 * mu_s[i] * (spin[j] * field[j] + @@ -472,7 +471,7 @@ double compute_demag_energy(fft_demag_plan *restrict plan, } -void finalize_plan(fft_demag_plan *restrict plan) { +void finalize_plan(fft_demag_plan * plan) { fftw_destroy_plan(plan->m_plan); fftw_destroy_plan(plan->h_plan); diff --git a/fidimag/common/dipolar/demag_oommf.c b/native/src/c_demag_oommf.cpp similarity index 92% rename from fidimag/common/dipolar/demag_oommf.c rename to native/src/c_demag_oommf.cpp index 8675802d..633b8184 100644 --- a/fidimag/common/dipolar/demag_oommf.c +++ b/native/src/c_demag_oommf.cpp @@ -1,9 +1,9 @@ -/* This file demag_oommf.c is taken from the OOMMF project (oommf/app/oxs/ext/demagcoef.cc +/* This file demag_oommf.c is taken from the OOMMF project (oommf/app/oxs/ext/demagcoef.cc downloaded from http://math.nist.gov/oommf/dist/oommf12a5rc_20120928.tar.gz) -with slightly modifications (change OC_REALWIDE to double) +with slightly modifications (change OC_REALWIDE to double) and distributed under this license shown below. */ -/* License +/* License OOMMF - Object Oriented MicroMagnetic Framework @@ -42,8 +42,48 @@ to copyright protection within the United States. #include #include -#include "dipolar.h" -#include "demagcoef.h" +#include "c_dipolar.h" +#include "c_demagcoef.h" +#include "c_constants.h" + +int AS_Compare(const void* px,const void* py) { + // Comparison based on absolute values + double x=fabs(*((const double *)px)); + double y=fabs(*((const double *)py)); + if(xy) return -1; + return 0; +} + + +double AccurateSum(int n,double *arr) { + // Order by decreasing magnitude + qsort(arr,n,sizeof(double), AS_Compare); + + // Add up using doubly compensated summation. If necessary, mark + // variables these "volatile" to protect against problems arising + // from extra precision. Also, don't expect the compiler to respect + // order with respect to parentheses at high levels of optimization, + // i.e., write "u=x; u-=(y-corr)" as opposed to "u=x-(y-corr)". + + double sum,corr,y,u,t,v,z,x,tmp; + + sum=arr[0]; corr=0; + for(int i=1;i 0. double result_sign=1.0; - if(x<0.0) - result_sign *= -1.0; - if(y<0.0) + if(x<0.0) + result_sign *= -1.0; + if(y<0.0) result_sign *= -1.0; x=fabs(x); y=fabs(y); z=fabs(z); // This function is even in z and /// odd in x and y. The fabs()'s simplify special case handling. double xsq=x*x,ysq=y*y,zsq=z*z; double R=xsq+ysq+zsq; - if(R<=0.0) + if(R<=0.0) return 0.0; else R=sqrt(R); @@ -369,7 +409,7 @@ double CalculateSDA01(double x,double y,double z, arr[26] = 8*Newell_g(x,y,z); - return AccurateSum(27,arr)/(4*WIDE_PI*l*h*e); + return AccurateSum(27,arr)/(4*M_PIl*l*h*e); /// factor Nxy from Newell's paper. } @@ -424,7 +464,7 @@ double DemagNxxAsymptotic(double x, double y, double z, term5 *= 0.25; } - double term7 = 0.0; + double term7 = 0.0; { double b1 = 32*dx4 - 40*dx2dy2 - 40*dx2dz2 + 12*dy4 + 10*dy2dz2 + 12*dz4; double b2 = -240*dx4 + 580*dx2dy2 + 20*dx2dz2 - 202*dy4 - 75*dy2dz2 + 22*dz4; @@ -442,7 +482,7 @@ double DemagNxxAsymptotic(double x, double y, double z, term7 *= (1.0/16.0); } - double Nxx = (-dx*dy*dz/(4*WIDE_PI)) * (((term7/R2 + term5)/R2 + term3)/(R2*R)); + double Nxx = (-dx*dy*dz/(4*M_PIl)) * (((term7/R2 + term5)/R2 + term3)/(R2*R)); // Error should be of order 1/R^9 return Nxx; @@ -505,7 +545,7 @@ double DemagNxyAsymptotic(double x, double y, double z, term7 *= (7.0/16.0); } - double Nxy = (-dx*dy*dz*x*y/(4*WIDE_PI*R2)) * (((term7/R2 + term5)/R2 + term3)/(R2*R)); + double Nxy = (-dx*dy*dz*x*y/(4*M_PIl*R2)) * (((term7/R2 + term5)/R2 + term3)/(R2*R)); // Error should be of order 1/R^9 return Nxy; diff --git a/native/src/c_energy.cpp b/native/src/c_energy.cpp new file mode 100644 index 00000000..5f5909b7 --- /dev/null +++ b/native/src/c_energy.cpp @@ -0,0 +1,34 @@ +#include "c_energy.h" +#include "c_constants.h" + +double Energy::compute_energy() { + double sum = 0; + for(int i = 0; i < n; i++) { + sum += energy[i]; + } + return sum * (dx * dy * dz * std::pow(unit_length, 3)); +} + +void Energy::_setup(MicroSim * sim) { + + this->nx = sim->nx; + this->ny = sim->ny; + this->nz = sim->nz; + this->n = sim->nx * ny * nz; + this->dx = sim->dx; + this->dy = sim->dy; + this->dz = sim->dz; + this->unit_length = sim->unit_length; + + // Arrays + this->spin = sim->spin; + this->Ms = sim->Ms; + this->Ms_inv = sim->Ms_inv; + this->coordinates = sim->coordinates; + this->ngbs = sim->ngbs; + this->energy = sim->energy; + this->field = sim->field; + // this->pins = sim->pins; + + set_up = true; +} diff --git a/fidimag/common/lib/llg.c b/native/src/c_llg.cpp similarity index 96% rename from fidimag/common/lib/llg.c rename to native/src/c_llg.cpp index ce438fdf..c7a743ee 100644 --- a/fidimag/common/lib/llg.c +++ b/native/src/c_llg.cpp @@ -1,4 +1,5 @@ -#include "common_clib.h" +#include "c_clib.h" +#include "c_vectormath.h" /* The right hand side of the LLG equation for the CVOde solver. This can be * used both for the micromagnetic and atomistic codes since m or S are unitless @@ -79,7 +80,7 @@ * value for the prefactor (6). */ -void llg_rhs(double *restrict dm_dt, double *restrict m, double *restrict h, double *restrict alpha, int *restrict pins, +void llg_rhs(double * dm_dt, double * m, double * h, double * alpha, int * pins, double gamma, int n, int do_precession, double default_c) { int i, j, k; @@ -166,7 +167,7 @@ void llg_rhs(double *restrict dm_dt, double *restrict m, double *restrict h, dou } -void llg_rhs_jtimes(double *restrict jtn, double *restrict m, double *restrict h, double *restrict mp, double *restrict hp, double *restrict alpha, int *restrict pins, +void llg_rhs_jtimes(double * jtn, double * m, double * h, double * mp, double * hp, double * alpha, int * pins, double gamma, int n, int do_precession, double default_c) { //#pragma omp parallel for private(i,j,k,coeff,mm, mh, c, hpi,hpj,hpk) diff --git a/native/src/c_micro_sim.cpp b/native/src/c_micro_sim.cpp new file mode 100644 index 00000000..37da3a33 --- /dev/null +++ b/native/src/c_micro_sim.cpp @@ -0,0 +1,47 @@ +#include "c_energy.h" +#include "c_micro_sim.h" +#include "c_constants.h" + + +void MicroSim::setup(int nx, int ny, int nz, double dx, double dy, double dz, + double unit_length, double *coordinates, int *ngbs, + double *spin, double *Ms, double *Ms_inv, + double *energy, double *field, int *pins + ) { + + this->nx = nx; + this->ny = ny; + this->nz = nz; + this->n = nx * ny * nz; + this->dx = dx; + this->dy = dy; + this->dz = dz; + this->unit_length = unit_length; + this->coordinates = coordinates; + this->ngbs = ngbs; + + this->spin = spin; + this->Ms = Ms; + this->Ms_inv = Ms_inv; + this->energy = energy; + this->field = field; + this->pins = pins; + + set_up = true; +} + +void MicroSim::add_interaction(Energy * interaction_ptr) { + + interactions.push_back(interaction_ptr); + // interactions_id.push_back(int_id); + +} + +void MicroSim::compute_effective_field(double t) { + + // Using indices + std::vector::iterator it; + for(it = interactions.begin(); it != interactions.end(); it++) { + (*it)->compute_field(t); + } +} diff --git a/fidimag/common/neb_method/nebm_geodesic_lib.c b/native/src/c_nebm_geodesic_lib.cpp similarity index 91% rename from fidimag/common/neb_method/nebm_geodesic_lib.c rename to native/src/c_nebm_geodesic_lib.cpp index 7bc28736..520245d9 100644 --- a/fidimag/common/neb_method/nebm_geodesic_lib.c +++ b/native/src/c_nebm_geodesic_lib.cpp @@ -1,10 +1,12 @@ -#include "nebm_geodesic_lib.h" -#include "nebm_lib.h" +#include "c_nebm_geodesic_lib.h" +#include "c_nebm_lib.h" #include "math.h" +#include "c_vectormath.h" -double compute_geodesic_Vincenty(double *restrict A, double *restrict B, + +double compute_geodesic_Vincenty(double * A, double * B, int n_dofs_image, - int *restrict material, + int * material, int n_dofs_image_material ) { @@ -68,9 +70,9 @@ double compute_geodesic_Vincenty(double *restrict A, double *restrict B, if (material[3 * i] > 0) { spin_i = 3 * i; - cross_product(A_cross_B, &A[spin_i], &B[spin_i]); + cross(A_cross_B, &A[spin_i], &B[spin_i]); A_cross_B_norm = compute_norm(A_cross_B, 3); - A_dot_B = dot_product(&A[spin_i], &B[spin_i], 3); + A_dot_B = dot(&A[spin_i], &B[spin_i], 3); geo_dist = atan2(A_cross_B_norm, A_dot_B); distance += geo_dist * geo_dist; @@ -82,9 +84,9 @@ double compute_geodesic_Vincenty(double *restrict A, double *restrict B, return distance; } -double compute_geodesic_GreatCircle(double *restrict A, double *restrict B, +double compute_geodesic_GreatCircle(double * A, double * B, int n_dofs_image, - int *restrict material, + int * material, int n_dofs_image_material ) { @@ -131,7 +133,7 @@ double compute_geodesic_GreatCircle(double *restrict A, double *restrict B, if (material[3 * i] > 0) { spin_i = 3 * i; - A_dot_B = dot_product(&A[spin_i], &B[spin_i], 3); + A_dot_B = dot(&A[spin_i], &B[spin_i], 3); // Prevent nans A_dot_B = fmax(-1.0, fmin(1.0, A_dot_B)); geo_dist = acos(A_dot_B); diff --git a/fidimag/common/neb_method/nebm_integrators.c b/native/src/c_nebm_integrators.cpp similarity index 83% rename from fidimag/common/neb_method/nebm_integrators.c rename to native/src/c_nebm_integrators.cpp index efaa39d3..79b2de34 100644 --- a/fidimag/common/neb_method/nebm_integrators.c +++ b/native/src/c_nebm_integrators.cpp @@ -1,13 +1,13 @@ #include "math.h" -#include "nebm_lib.h" +#include "c_nebm_lib.h" #include #include -double step_Verlet_C(double *restrict forces, - double *restrict forces_prev, - double *restrict velocities, - double *restrict velocities_new, - double *restrict y, +double step_Verlet_C(double * forces, + double * forces_prev, + double * velocities, + double * velocities_new, + double * y, double t, double h, double mass, @@ -17,8 +17,6 @@ double step_Verlet_C(double *restrict forces, ) { - printf("Hellooooooo\n"); - update_field(t, y); int im_idx; diff --git a/fidimag/common/neb_method/nebm_lib.c b/native/src/c_nebm_lib.cpp similarity index 82% rename from fidimag/common/neb_method/nebm_lib.c rename to native/src/c_nebm_lib.cpp index 8bb0611c..6ac761b3 100644 --- a/fidimag/common/neb_method/nebm_lib.c +++ b/native/src/c_nebm_lib.cpp @@ -1,71 +1,10 @@ -#include "nebm_lib.h" +#include "c_nebm_lib.h" #include "math.h" #include +#include "c_vectormath.h" -void cross_product(double *restrict output, double * A, double * B){ - /* Cross product between arrays A and B, assuming they - * have length = 3 - * Every resulting component is stored in the *output array - */ - - output[0] = A[1] * B[2] - A[2] * B[1]; - output[1] = A[2] * B[0] - A[0] * B[2]; - output[2] = A[0] * B[1] - A[1] * B[0]; -} -double dot_product(double * A, double * B, int n){ - /* Dot product between arrays A and B , assuming they - * have length n */ - - double dotp = 0; - for(int i = 0; i < n; i++){ - dotp += A[i] * B[i]; - } - - return dotp; -} - -double compute_norm(double *restrict a, int n) { - /* Compute the norm of an array *a - * - * ARGUMENTS: - - * n :: length of a - - */ - - double norm = 0; - for(int i = 0; i < n; i++){ norm += a[i] * a[i]; } - norm = sqrt(norm); - return norm; -} - -void normalise(double *restrict a, int n){ - - /* Normalise the *a array, whose length is n (3 * number of nodes in - * cartesian, and 2 * number of nodes in spherical) To do this we compute - * the length of *a : - * - * SQRT[ a[0] ** 2 + a[1] ** 2 + ... ] - * - * and divide every *a entry by that length - */ - - double length; - - length = compute_norm(a, n); - - if (length > 0){ - length = 1.0 / length; - } - - for(int i = 0; i < n; i++){ - a[i] *= length; - } -} - - -void normalise_images_C(double *restrict y, int n_images, int n_dofs_image){ +void normalise_images_C(double * y, int n_images, int n_dofs_image){ int i; @@ -80,7 +19,8 @@ void normalise_images_C(double *restrict y, int n_images, int n_dofs_image){ } } -void normalise_spins_C(double *restrict y, int n_images, int n_dofs_image){ + +void normalise_spins_C(double * y, int n_images, int n_dofs_image){ int i, j; @@ -102,7 +42,7 @@ void normalise_spins_C(double *restrict y, int n_images, int n_dofs_image){ /* ------------------------------------------------------------------------- */ -void compute_tangents_C(double *restrict tangents, double *restrict y, double *restrict energies, +void compute_tangents_C(double * tangents, double * y, double * energies, int n_dofs_image, int n_images ) { @@ -167,8 +107,8 @@ void compute_tangents_C(double *restrict tangents, double *restrict y, double *r double *t_plus; double *t_minus; - t_plus = malloc(n_dofs_image * sizeof(double)); - t_minus = malloc(n_dofs_image * sizeof(double)); + t_plus = new double[n_dofs_image]; + t_minus = new double[n_dofs_image]; double deltaE_plus, deltaE_minus; double deltaE_MAX, deltaE_MIN; @@ -252,21 +192,21 @@ void compute_tangents_C(double *restrict tangents, double *restrict y, double *r /* ----------------------------------------------------------------- */ } // Close loop in images - free(t_plus); - free(t_minus); + delete[] t_plus; + delete[] t_minus; } // Close main function /* ------------------------------------------------------------------------- */ void compute_spring_force_C( - double *restrict spring_force, - double *restrict y, - double *restrict tangents, - double *restrict k, + double * spring_force, + double * y, + double * tangents, + double * k, int n_images, int n_dofs_image, - double *restrict distances + double * distances ) { /* Compute the spring force for every image of an energy band, which is @@ -321,11 +261,11 @@ void compute_spring_force_C( /* ------------------------------------------------------------------------- */ -void compute_effective_force_C(double *restrict G, - double *restrict tangents, - double *restrict gradientE, - double *restrict spring_force, - int *restrict climbing_image, +void compute_effective_force_C(double * G, + double * tangents, + double * gradientE, + double * spring_force, + int * climbing_image, int n_images, int n_dofs_image) { @@ -365,7 +305,7 @@ void compute_effective_force_C(double *restrict G, double * gradE = &gradientE[im_idx]; double * sf = &spring_force[im_idx]; - gradE_dot_t = dot_product(gradE, t, n_dofs_image); + gradE_dot_t = dot(gradE, t, n_dofs_image); if(climbing_image[i] == 0) { for(j = 0; j < n_dofs_image; j++) { @@ -395,9 +335,9 @@ void compute_effective_force_C(double *restrict G, // arrays. The *distances array saves the distances as: // [ |Y1 - Y0|, |Y2 - Y1|, ... ] // Path distances are the total distances relative to the 0-th image -void compute_image_distances(double *restrict distances, - double *restrict path_distances, - double *restrict y, +void compute_image_distances(double * distances, + double * path_distances, + double * y, int n_images, int n_dofs_image, double (* compute_distance)(double *, @@ -405,7 +345,7 @@ void compute_image_distances(double *restrict distances, int, int *, int), - int *restrict material, + int * material, int n_dofs_image_material ) { @@ -427,7 +367,7 @@ void compute_image_distances(double *restrict distances, // ---------------------------------------------------------------------------- -void project_vector_C(double *restrict vector, double *restrict y_i, +void project_vector_C(double * vector, double * y_i, int n_dofs_image ){ @@ -464,13 +404,13 @@ void project_vector_C(double *restrict vector, double *restrict y_i, // dot product of the i-th vector (i.e. using // ( vector[j], vector[j+1], vector[j+2] ) , j=0,3,6, ... ) // and then compute the projection for the 3 components of the vector - if (j % 3 == 0) v_dot_m_i = dot_product(&vector[j], &y_i[j], 3); + if (j % 3 == 0) v_dot_m_i = dot(&vector[j], &y_i[j], 3); vector[j] = vector[j] - v_dot_m_i * y_i[j]; } } -void project_images_C(double *restrict vector, double *restrict y, +void project_images_C(double * vector, double * y, int n_images, int n_dofs_image ){ @@ -520,8 +460,8 @@ void project_images_C(double *restrict vector, double *restrict y, // ---------------------------------------------------------------------------- -double compute_distance_cartesian(double *restrict A, double *restrict B, int n_dofs_image, - int *restrict material, int n_dofs_image_material +double compute_distance_cartesian(double * A, double * B, int n_dofs_image, + int * material, int n_dofs_image_material ) { /* Compute the distance between two images, A and B, discarding the sites @@ -568,8 +508,8 @@ double compute_distance_cartesian(double *restrict A, double *restrict B, int n_ * */ -void compute_dYdt(double *restrict Y, double *restrict G, double *restrict dYdt, - int *restrict pins, +void compute_dYdt(double * Y, double * G, double * dYdt, + int * pins, int n_dofs_image ) { @@ -584,8 +524,8 @@ void compute_dYdt(double *restrict Y, double *restrict G, double *restrict dYdt, continue; } - double Y_dot_Y = Y[j] * Y[j] + Y[j + 1] * Y[j + 1] + Y[j + 2] * Y[j + 2]; - double Y_dot_G = Y[j] * G[j] + Y[j + 1] * G[j + 1] + Y[j + 2] * G[j + 2]; + double Y_dot_Y = dot(&Y[j], &Y[j], 3); + double Y_dot_G = dot(&Y[j], &G[j], 3); // (Y * Y) G - (Y * G) Y = - Y x (Y x G) dYdt[j] = Y_dot_Y * G[j] - Y_dot_G * Y[j]; dYdt[j + 1] = Y_dot_Y * G[j + 1] - Y_dot_G * Y[j + 1]; @@ -602,7 +542,7 @@ void compute_dYdt(double *restrict Y, double *restrict G, double *restrict dYdt, } } -void compute_dYdt_C(double *restrict y, double *restrict G, double *restrict dYdt, int *restrict pins, +void compute_dYdt_C(double * y, double * G, double * dYdt, int * pins, int n_images, int n_dofs_image) { #pragma omp parallel for schedule(static) for(int i = 1; i < n_images - 1; i++){ @@ -617,8 +557,8 @@ void compute_dYdt_C(double *restrict y, double *restrict G, double *restrict dYd // dYdt functions without the magic correction factor // Necessary to implement more standard integrators, such as Verlet -void compute_dYdt_nc(double *restrict Y, double *restrict G, double *restrict dYdt, - int *restrict pins, +void compute_dYdt_nc(double * Y, double * G, double * dYdt, + int * pins, int n_dofs_image ) { @@ -633,17 +573,18 @@ void compute_dYdt_nc(double *restrict Y, double *restrict G, double *restrict dY continue; } - double Y_dot_Y = Y[j] * Y[j] + Y[j + 1] * Y[j + 1] + Y[j + 2] * Y[j + 2]; - double Y_dot_G = Y[j] * G[j] + Y[j + 1] * G[j + 1] + Y[j + 2] * G[j + 2]; + double Y_dot_Y = dot(&Y[j], &Y[j], 3); + double Y_dot_G = dot(&Y[j], &G[j], 3); + // (Y * Y) G - (Y * G) Y = - Y x (Y x G) - dYdt[j] = Y_dot_Y * G[j] - Y_dot_G * Y[j]; - dYdt[j + 1] = Y_dot_Y * G[j + 1] - Y_dot_G * Y[j + 1]; - dYdt[j + 2] = Y_dot_Y * G[j + 2] - Y_dot_G * Y[j + 2]; + dYdt[j] = Y_dot_Y * G[j] - Y_dot_G * Y[j]; + dYdt[j + 1] = Y_dot_Y * G[j + 1] - Y_dot_G * Y[j + 1]; + dYdt[j + 2] = Y_dot_Y * G[j + 2] - Y_dot_G * Y[j + 2]; } } -void compute_dYdt_nc_C(double *restrict y, double *restrict G, double *restrict dYdt, - int *restrict pins, +void compute_dYdt_nc_C(double * y, double * G, double * dYdt, + int * pins, int n_images, int n_dofs_image) { #pragma omp parallel for schedule(static) for(int i = 1; i < n_images - 1; i++){ diff --git a/fidimag/common/neb_method/nebm_spherical_lib.c b/native/src/c_nebm_spherical_lib.cpp similarity index 78% rename from fidimag/common/neb_method/nebm_spherical_lib.c rename to native/src/c_nebm_spherical_lib.cpp index 7ef9641f..4a0d86c1 100644 --- a/fidimag/common/neb_method/nebm_spherical_lib.c +++ b/native/src/c_nebm_spherical_lib.cpp @@ -1,25 +1,26 @@ -#include "nebm_spherical_lib.h" -#include "nebm_lib.h" +#include "c_nebm_spherical_lib.h" +#include "c_nebm_lib.h" #include "math.h" +#include "c_constants.h" -double compute_norm_spherical(double *restrict a, int n, int scale) { +double compute_norm_spherical(double * a, int n, int scale) { /* Compute the norm of an array *a. *a is assumed to have spherical * coordinates as: - + * a = [ theta_0 phi_0 theta_1 phi_1 ... ] - + * This function is necessary to either normalise a vector or compute the * distance between two images (see below) * * To compute the norm, we redefine the angles, so they always lie on the - * [-PI, PI] range, rather than [0, 2PI] for phi. + * [-PI, PI] range, rather than [0, 2PI] for phi. * * ARGUMENTS: - + * n :: length of a * * scale :: If scale is zero, we do not scale the norm - + * Notice that, when computing the DISTANCE between two images, Y_i, Y_i+1 * for example, we just do the difference (Y_i - Y_(i+1)) and then compute * its norm (using this function) but RESCALING the length by the number @@ -32,10 +33,10 @@ double compute_norm_spherical(double *restrict a, int n, int scale) { for(int i = 0; i < n; i++){ - if (a[i] > WIDE_PI){ - a[i] = 2 * WIDE_PI - a[i]; - } else if(a[i] < -WIDE_PI){ - a[i] += 2 * WIDE_PI; + if (a[i] > M_PIl){ + a[i] = 2 * M_PIl - a[i]; + } else if(a[i] < -M_PIl){ + a[i] += 2 * M_PIl; } norm += a[i] * a[i]; @@ -47,7 +48,7 @@ double compute_norm_spherical(double *restrict a, int n, int scale) { return norm; } -void normalise_spherical(double *restrict a, int n){ +void normalise_spherical(double * a, int n){ /* Normalise the *a array, whose length is n (3 * number of nodes in * cartesian, and 2 * number of nodes in spherical) To do this we compute @@ -72,7 +73,7 @@ void normalise_spherical(double *restrict a, int n){ } -void normalise_images_spherical_C(double *restrict y, int n_images, int n_dofs_image){ +void normalise_images_spherical_C(double * y, int n_images, int n_dofs_image){ int i; @@ -89,8 +90,8 @@ void normalise_images_spherical_C(double *restrict y, int n_images, int n_dofs_i /* ------------------------------------------------------------------------- */ -double compute_distance_spherical(double *restrict A, double *restrict B, int n, - int *restrict material, int n_dofs_image_material +double compute_distance_spherical(double * A, double * B, int n, + int * material, int n_dofs_image_material ) { double A_minus_B[n_dofs_image_material]; diff --git a/fidimag/common/lib/steepest_descent.c b/native/src/c_steepest_descent.cpp similarity index 99% rename from fidimag/common/lib/steepest_descent.c rename to native/src/c_steepest_descent.cpp index f1d860b0..a7d19240 100644 --- a/fidimag/common/lib/steepest_descent.c +++ b/native/src/c_steepest_descent.cpp @@ -1,5 +1,6 @@ -#include "common_clib.h" +#include "c_clib.h" #include "math.h" +#include "c_vectormath.h" void sd_update_spin (double *spin, double *spin_last, double *magnetisation, double *mxH, double *mxmxH, double *mxmxH_last, double tau, diff --git a/fidimag/common/lib/stt.c b/native/src/c_stt.cpp similarity index 93% rename from fidimag/common/lib/stt.c rename to native/src/c_stt.cpp index b2faeb10..f76fa06f 100644 --- a/fidimag/common/lib/stt.c +++ b/native/src/c_stt.cpp @@ -1,4 +1,5 @@ -#include "common_clib.h" +#include "c_clib.h" +#include "c_vectormath.h" /* Calculation of the STT field using the Zhang Li formalism, which * is @@ -51,8 +52,8 @@ * Neighbouring sites where there is no material, has index -1 * */ -void compute_stt_field_c(double *restrict spin, double *restrict field, double *restrict jx, double *restrict jy, double *restrict jz, - double dx, double dy, double dz, int *restrict ngbs, int n) { +void compute_stt_field_c(double * spin, double * field, double * jx, double * jy, double * jz, + double dx, double dy, double dz, int * ngbs, int n) { //#pragma omp parallel for for (int i = 0; i < 3 * n; i++) { @@ -159,8 +160,8 @@ void compute_stt_field_c(double *restrict spin, double *restrict field, double * } -void llg_stt_rhs(double *restrict dm_dt, double *restrict m, double *restrict h, double *restrict h_stt, - double *restrict alpha, double beta, double u0, double gamma, int n) { +void llg_stt_rhs(double * dm_dt, double * m, double * h, double * h_stt, + double * alpha, double beta, double u0, double gamma, int n) { #pragma omp parallel for for (int index = 0; index < n; index++) { @@ -220,8 +221,8 @@ void llg_stt_rhs(double *restrict dm_dt, double *restrict m, double *restrict h, } -void llg_stt_cpp(double *restrict dm_dt, double *restrict m, double *restrict h, double *restrict p, - double *restrict alpha, int *restrict pins, double *restrict a_J, double beta, double gamma, int n) { +void llg_stt_cpp(double * dm_dt, double * m, double * h, double * p, + double * alpha, int * pins, double * a_J, double beta, double gamma, int n) { #pragma omp parallel for for (int index = 0; index < n; index++) { diff --git a/fidimag/common/dipolar/tensor_2dpbc.c b/native/src/c_tensor_2dpbc.cpp similarity index 93% rename from fidimag/common/dipolar/tensor_2dpbc.c rename to native/src/c_tensor_2dpbc.cpp index d9104908..db9acb9a 100644 --- a/fidimag/common/dipolar/tensor_2dpbc.c +++ b/native/src/c_tensor_2dpbc.cpp @@ -1,8 +1,9 @@ #include #include -#include "dipolar.h" -#include "demagcoef.h" -#include "tensor_2dpbc.h" +#include "c_dipolar.h" +#include "c_demagcoef.h" +#include "c_tensor_2dpbc.h" +#include "c_constants.h" double DemagTensorNormal(enum Type_Nij comp,double x,double y,double z,double a,double b,double c) { @@ -29,13 +30,13 @@ double DemagTensorAsymptotic(enum Type_Nij comp,double x,double y,double z,doubl case Tensor_xz:return DemagNxyAsymptotic(x,z,y,a,c,b); case Tensor_yz:return DemagNxyAsymptotic(y,z,x,b,c,a); } - + return -1; } double DemagTensorDipolar(enum Type_Nij comp,double x,double y,double z){ - + switch(comp){ case Tensor_xx:return DemagNxxDipolar(x,y,z); case Tensor_yy:return DemagNxxDipolar(y,x,z); @@ -66,12 +67,12 @@ double DemagTensorInfinite(enum Type_Nij comp, double x,double y,double z,double double compute_single_tensor(enum Type_Nij comp, int g, double x, double y, double z, - double a, double b, double c, int xdim, int ydim, + double a, double b, double c, int xdim, int ydim, double asymptotic_radius_sq, double dipolar_radius_sq){ if ((comp == Tensor_xy || comp == Tensor_xz || comp == Tensor_yz) && x * y == 0) return 0.0; - double Tx = xdim*a, Ty = ydim*b, cof2 = a * b * c / (4 * WIDE_PI); + double Tx = xdim*a, Ty = ydim*b, cof2 = a * b * c / (4 * M_PIl); double* tmpx = (double *)malloc(sizeof(double)*(2 * g + 2)); double* tmpy = (double *)malloc(sizeof(double)*(2 * g + 2)); @@ -107,12 +108,12 @@ double compute_single_tensor(enum Type_Nij comp, int g, double x, double y, doub } double compute_single_tensor_finitely(enum Type_Nij comp, int gx, int gy, double x, double y, double z, - double a, double b, double c, int xdim, int ydim, + double a, double b, double c, int xdim, int ydim, double asymptotic_radius_sq, double dipolar_radius_sq){ if ((comp == Tensor_xy || comp == Tensor_xz || comp == Tensor_yz) && x * y == 0) return 0.0; - double Tx = xdim*a, Ty = ydim*b, cof2 = a * b * c / (4 * WIDE_PI); + double Tx = xdim*a, Ty = ydim*b, cof2 = a * b * c / (4 * M_PIl); double* tmpx = (double *)malloc(sizeof(double)*(2 * gx + 1)); double* tmpy = (double *)malloc(sizeof(double)*(2 * gy + 1)); double tpx, tpy, radius_sq; @@ -152,20 +153,20 @@ int FindG(enum Type_Nij comp, double v, double Tx, double Ty, double pbc_2d_erro case Tensor_xz: case Tensor_yz: case Tensor_xx: - tmp = v / (4 * WIDE_PI * (Tx * Tx) * sqrt(Tx * Tx + Ty * Ty) * pbc_2d_error); + tmp = v / (4 * M_PIl * (Tx * Tx) * sqrt(Tx * Tx + Ty * Ty) * pbc_2d_error); return (int) pow(tmp, 1 / 3.0) + 1; case Tensor_yy: - tmp = v / (4 * WIDE_PI * (Ty * Ty) * sqrt(Tx * Tx + Ty * Ty) * pbc_2d_error); + tmp = v / (4 * M_PIl * (Ty * Ty) * sqrt(Tx * Tx + Ty * Ty) * pbc_2d_error); return (int) pow(tmp, 1 / 3.0) + 1; case Tensor_zz: - tmp = v * sqrt(Tx * Tx + Ty * Ty) / (4 * WIDE_PI * (Tx * Ty * Tx * Ty) * pbc_2d_error); + tmp = v * sqrt(Tx * Tx + Ty * Ty) / (4 * M_PIl * (Tx * Ty * Tx * Ty) * pbc_2d_error); return (int) pow(tmp, 1 / 3.0) + 1; } return 0; } void fill_demag_tensors_c(fft_demag_plan *plan, double *tensors){ - + int nx = plan->nx; int ny = plan->ny; @@ -179,7 +180,7 @@ void fill_demag_tensors_c(fft_demag_plan *plan, double *tensors){ int len_xy = lenx * leny; int x,y,z, id, index; - + for (int k = 0; k < lenz; k++) { for (int j = 0; j < leny; j++) { for (int i = 0; i < lenx; i++) { @@ -209,7 +210,7 @@ void fill_demag_tensors_c(fft_demag_plan *plan, double *tensors){ //compute the demag tensors, i.e, H=-N.M void compute_demag_tensors_2dpbc(fft_demag_plan *plan, double *tensors, double pbc_2d_error, int sample_repeat_nx, int sample_repeat_ny, double dipolar_radius){ - + int gxx = sample_repeat_nx; int gyy = sample_repeat_ny; int gzz = 0; @@ -266,7 +267,7 @@ void compute_demag_tensors_2dpbc(fft_demag_plan *plan, double *tensors, gyy = FindG(Tensor_yy, dx * dy*dz, nx*dx, ny * dy, pbc_2d_error); gzz = FindG(Tensor_zz, dx * dy*dz, nx*dx, ny * dy, pbc_2d_error); //printf("gxx=%d gyy=%d gzz=%d\n",gxx,gyy,gzz); - + for (k = 0; k < nz; k++) { for (j = 0; j < ny; j++) { for (i = 0; i < nx; i++) { @@ -282,17 +283,12 @@ void compute_demag_tensors_2dpbc(fft_demag_plan *plan, double *tensors, tensors[id+3*nxyz] = compute_single_tensor(Tensor_xy, gxx, x, y, z, dx, dy, dz, nx, ny, asymptotic_radius_sq, dipolar_radius_sq); tensors[id+4*nxyz] = compute_single_tensor(Tensor_xz, gxx, x, y, z, dx, dy, dz, nx, ny, asymptotic_radius_sq, dipolar_radius_sq); tensors[id+5*nxyz] = compute_single_tensor(Tensor_yz, gxx, x, y, z, dx, dy, dz, nx, ny, asymptotic_radius_sq, dipolar_radius_sq); - + //printf("x=%g y=%g z=%g nxx=%g nyy=%g nzz=%g\n",x,y,z,plan->tensor_xx[id],plan->tensor_yy[id],plan->tensor_zz[id]); } } } } - -} - - - - +} diff --git a/native/src/c_vectormath.cpp b/native/src/c_vectormath.cpp new file mode 100644 index 00000000..4737a44c --- /dev/null +++ b/native/src/c_vectormath.cpp @@ -0,0 +1,102 @@ +#include +#include "c_vectormath.h" + + +double cross_x(double a0, double a1, double a2, double b0, double b1, double b2) { + return a1*b2 - a2*b1; +} + +double cross_y(double a0, double a1, double a2, double b0, double b1, double b2) { + return a2*b0 - a0*b2; +} + +double cross_z(double a0, double a1, double a2, double b0, double b1, double b2) { + return a0*b1 - a1*b0; +} + +void normalise(double *a, int n) { + /* Normalise the *a array, whose length is n (3 * number of nodes in + * cartesian, and 2 * number of nodes in spherical) To do this we compute + * the length of *a : + * + * SQRT[ a[0] ** 2 + a[1] ** 2 + ... ] + * + * and divide every *a entry by that length + */ + + double length; + + length = compute_norm(a, n); + + if (length > 0){ + length = 1.0 / length; + } + + for(int i = 0; i < n; i++){ + a[i] *= length; + } +} + +void normalise(double *m, int *pins, int n) { + // Normalises n three vectors + int i, j, k; + double mdotm; + for (int id = 0; id < n; id++) { + i = 3*id; + j = i + 1; + k = j + 1; + + if (pins[id]>0) continue; + + mdotm = std::sqrt(m[i] * m[i] + m[j] * m[j] + m[k] * m[k]); + if(mdotm > 0) { + mdotm = 1 / mdotm; + m[i] *= mdotm; + m[j] *= mdotm; + m[k] *= mdotm; + } + } +} + + + + +void cross(double * output, double * A, double * B){ + /* Cross product between arrays A and B, assuming they + * have length = 3 + * Every resulting component is stored in the *output array + */ + + output[0] = A[1] * B[2] - A[2] * B[1]; + output[1] = A[2] * B[0] - A[0] * B[2]; + output[2] = A[0] * B[1] - A[1] * B[0]; +} + +double dot(double * A, double * B, int n){ + /* Dot product between arrays A and B , assuming they + * have length n */ + + double dotp = 0; + for(int i = 0; i < n; i++){ + dotp += A[i] * B[i]; + } + + return dotp; +} + +double compute_norm(double * a, int n) { + /* Compute the norm of an array *a + * + * ARGUMENTS: + + * n :: length of a + + */ + + double norm = 0; + for(int i = 0; i < n; i++) { + norm += a[i] * a[i]; + } + norm = std::sqrt(norm); + return norm; +} diff --git a/native/src/m_anis.cpp b/native/src/m_anis.cpp new file mode 100644 index 00000000..0bed2b37 --- /dev/null +++ b/native/src/m_anis.cpp @@ -0,0 +1,58 @@ +#include "m_clib.h" +#include "c_energy.h" +#include "c_constants.h" + +// TO BE DEPRECATED:: + +void compute_uniaxial_anis(double * m, double * field, double * energy, double * Ms_inv, + double * Ku, double * axis, int nx, int ny, int nz) { + + int n = nx * ny * nz; + +#pragma omp parallel for + for (int i = 0; i < n; i++) { + int j = 3 * i; + + if (Ms_inv[i] == 0.0){ + field[j] = 0; + field[j + 1] = 0; + field[j + 2] = 0; + energy[i] = 0; + continue; + } + + double m_u = m[j] * axis[j] + m[j + 1] * axis[j + 1] + m[j + 2] * axis[j + 2]; + + field[j] = 2 * Ku[i] * m_u * Ms_inv[i] * MU0_INV * axis[j]; + field[j + 1] = 2 * Ku[i] * m_u * Ms_inv[i] * MU0_INV * axis[j + 1]; + field[j + 2] = 2 * Ku[i] * m_u * Ms_inv[i] * MU0_INV * axis[j + 2]; + + energy[i] = Ku[i] * (1 - m_u * m_u); + } +} + +// ---------------------------------------------------------------------------- + +void AnisotropyEnergy::compute_field(double t) { + + for (int i = 0; i < n; i++) { + + int j = 3 * i; + + if (Ms_inv[i] == 0.0){ + field[j] = 0; + field[j + 1] = 0; + field[j + 2] = 0; + energy[i] = 0; + continue; + } + + double m_u = spin[j] * axis[j] + spin[j + 1] * axis[j + 1] + spin[j + 2] * axis[j + 2]; + + field[j] = 2 * Ku[i] * m_u * Ms_inv[i] * MU0_INV * axis[j]; + field[j + 1] = 2 * Ku[i] * m_u * Ms_inv[i] * MU0_INV * axis[j + 1]; + field[j + 2] = 2 * Ku[i] * m_u * Ms_inv[i] * MU0_INV * axis[j + 2]; + + energy[i] = Ku[i] * (1 - m_u * m_u); + } +} diff --git a/fidimag/micro/lib/baryakhtar/llg.c b/native/src/m_baryakhtar_llg.cpp similarity index 96% rename from fidimag/micro/lib/baryakhtar/llg.c rename to native/src/m_baryakhtar_llg.cpp index b1f0b0cf..04c1356f 100644 --- a/fidimag/micro/lib/baryakhtar/llg.c +++ b/native/src/m_baryakhtar_llg.cpp @@ -1,4 +1,7 @@ -#include "baryakhtar_clib.h" +#include "c_vectormath.h" +#include "m_clib.h" +#include "m_baryakhtar_clib.h" + void llg_rhs_baryakhtar(double *dm_dt, double *m, double *h, double *delta_h, double *alpha, double beta, int *pins, diff --git a/fidimag/micro/lib/baryakhtar/util.c b/native/src/m_baryakhtar_util.cpp similarity index 98% rename from fidimag/micro/lib/baryakhtar/util.c rename to native/src/m_baryakhtar_util.cpp index a0e82da2..4b4a9f60 100644 --- a/fidimag/micro/lib/baryakhtar/util.c +++ b/native/src/m_baryakhtar_util.cpp @@ -1,4 +1,4 @@ -#include "baryakhtar_clib.h" +#include "m_baryakhtar_clib.h" void compute_laplace_m(double *m, double *field, double *Ms, double dx, double dy, double dz, int nx, int ny, int nz) { diff --git a/fidimag/micro/lib/dmi.c b/native/src/m_dmi.cpp similarity index 75% rename from fidimag/micro/lib/dmi.c rename to native/src/m_dmi.cpp index bdb1e7b3..ec0dd209 100644 --- a/fidimag/micro/lib/dmi.c +++ b/native/src/m_dmi.cpp @@ -1,7 +1,10 @@ -#include "micro_clib.h" +#include "c_vectormath.h" +#include "m_clib.h" #include #include #include +#include "c_energy.h" +#include "c_constants.h" /* Functions to Compute the micromagnetic Dzyaloshinskii Moriya interaction * field and energy using the @@ -266,12 +269,12 @@ * the last 6 terms are set to zero * */ -void dmi_field(double *restrict m, double *restrict field, - double *restrict energy, double *restrict Ms_inv, - double *restrict D, int n_dmis, +void dmi_field(double * m, double * field, + double * energy, double * Ms_inv, + double * D, int n_dmis, double *dmi_vector, double dx, double dy, double dz, - int n, int *restrict ngbs) { + int n, int * ngbs) { /* These are for the DMI prefactor or coefficient */ double dxs[6] = {dx, dx, dy, dy, dz, dz}; @@ -363,3 +366,101 @@ void dmi_field(double *restrict m, double *restrict field, field[3 * i + 2] = fz * Ms_inv[i] * MU0_INV; } } + +// ---------------------------------------------------------------------------- +// Using the new Energy class:: + + +void DMIEnergy::compute_field(double t) { + + // Should use: std::unique_ptr dxs ?? + /* These are for the DMI prefactor or coefficient */ + double dxs[6] = {dx, dx, dy, dy, dz, dz}; + + /* Here we iterate through every mesh node */ + for (int i = 0; i < n; i++) { + double DMIc; + double fx = 0, fy = 0, fz = 0; + int idnm; // Index for the magnetisation matrix + int idn = 6 * i; // index for the neighbours + /* Set a zero field for sites without magnetic material */ + if (Ms_inv[i] == 0.0){ + field[3 * i] = 0; + field[3 * i + 1] = 0; + field[3 * i + 2] = 0; + energy[i] = 0; + continue; + } + + /* Here we iterate through the neighbours. Remember: + * j = 0, 1, 2, 3, 4, 5 --> -x, +x, -y, +y, -z, +z + Previously we had n_dmi_neighbours as a parameter, + but now we just skip if the vector entries are zero + instead, because we can have cases such as in D_n + where the D1 component has -x,+x,-y,+y,0,0 + and the D2 component has 0,0,0,0,-z,+z. + */ + for (int j = 0; j < 6; j++) { + /* Add the DMI field x times for every DMI constant */ + for (int k = 0; k < n_dmis; k++) { + + // starting index of the DMI vector for this neighbour (j) + // (remember we have 18 comps of dmi_vector per DMI constant) + int ngbr_idx_D = k * 18 + 3 * j; + + /* We skip the neighbour if + (a) it doesn't exist (ngbs[idn + j] = -1) + (b) there is no material there + (c) DMI value is zero there + */ + if ((ngbs[idn + j] != -1) && (Ms_inv[ngbs[idn + j]] != 0)) { + /* We do here: -(D / dx_i) * ( r_{ij} X M_{j} ) + * The cross_i function gives the i component of the + * cross product. The coefficient is computed according + * to the DMI strength of the current lattice site. + * For the denominator, for example, if j=2 or 3, then + * dxs[j] = dy + */ + + /* check the DMI vector exists */ + if (dmi_vector[ngbr_idx_D ] != 0 || + dmi_vector[ngbr_idx_D + 1] != 0 || + dmi_vector[ngbr_idx_D + 2] != 0 ) { + /* Notice the x component of the cross product of +-x + * times anything is zero (similar for the other comps) */ + + // Get DMI constant from neighbour and scale + // DMI components are in consecutive order per neighbour, + // i.e. if 2 DMI consts: [D0 D1 D0 D1 .... ] + DMIc = -D[n_dmis * ngbs[idn + j] + k] / dxs[j]; + + idnm = 3 * ngbs[idn + j]; // index for magnetisation + + fx += DMIc * cross_x(dmi_vector[ngbr_idx_D], + dmi_vector[ngbr_idx_D + 1], + dmi_vector[ngbr_idx_D + 2], + spin[idnm], spin[idnm + 1], spin[idnm + 2]); + fy += DMIc * cross_y(dmi_vector[ngbr_idx_D], + dmi_vector[ngbr_idx_D + 1], + dmi_vector[ngbr_idx_D + 2], + spin[idnm], spin[idnm + 1], spin[idnm + 2]); + fz += DMIc * cross_z(dmi_vector[ngbr_idx_D], + dmi_vector[ngbr_idx_D + 1], + dmi_vector[ngbr_idx_D + 2], + spin[idnm], spin[idnm + 1], spin[idnm + 2]); + + } + } + } // Close for loop through n of DMI constants + } // Close for loop through neighbours per mesh site + + /* Energy as: (-mu0 * Ms / 2) * [ H_dmi * m ] */ + energy[i] = -0.5 * (fx * spin[3 * i] + fy * spin[3 * i + 1] + + fz * spin[3 * i + 2]); + + /* Update the field H_dmi which has the same structure than *m */ + field[3 * i] = fx * Ms_inv[i] * MU0_INV; + field[3 * i + 1] = fy * Ms_inv[i] * MU0_INV; + field[3 * i + 2] = fz * Ms_inv[i] * MU0_INV; + } +} diff --git a/native/src/m_driver.cpp b/native/src/m_driver.cpp new file mode 100644 index 00000000..91229f0d --- /dev/null +++ b/native/src/m_driver.cpp @@ -0,0 +1,244 @@ +#include "m_driver.h" +#include "c_vectormath.h" +#include + + +// The RHS function to integrate the LLG equation +// Extra parameters are passed to the integrator in a struct +void LLG_RHS(double * m, std::vector& dmdt, unsigned int n, double t, + LLG_params * llg_params) { + + // Only if using void pointer for parameters: + // LLG_params * llg_params = static_cast(parameters); + + // Compute LLG RHS in m + unsigned int i, j, k; + double coeff, mm, mh, c; + double hpi, hpj, hpk; + + // #pragma omp parallel for private(i,j,k,coeff,mm, mh, c, hpi,hpj,hpk) + for (int id = 0; id < n; id++) { + // Indexes for the 3 components of the spin (magnetic moment) + // at the i-th lattice (mesh) site --> x, y, z + i = 3 * id; + j = i + 1; + k = i + 2; + + // Pinned spins do not follow the dynamical equation + if (llg_params->pins[id] > 0) { + dmdt[i] = 0; + dmdt[j] = 0; + dmdt[k] = 0; + continue; + } + + coeff = -llg_params->gamma / (1.0 + llg_params->alpha[id] * llg_params->alpha[id]); + + // Dot products + mm = m[i] * m[i] + m[j] * m[j] + m[k] * m[k]; + mh = m[i] * llg_params->field[i] + + m[j] * llg_params->field[j] + + m[k] * llg_params->field[k]; + + // Usually, m is normalised, i.e., mm=1; + // so hp = mm.h - mh.m = -m x (m x h) + // We set here the perpendicular componenet of the field + // but using the (m * m) product + hpi = mm * llg_params->field[i] - mh * m[i]; + hpj = mm * llg_params->field[j] - mh * m[j]; + hpk = mm * llg_params->field[k] - mh * m[k]; + + // IMPORTAthis->NT: do not ignore mm !! + // What we've found is that if we igonre mm, i.e. using + // hpi = h[i] - mh * m[i]; + // hpj = h[j] - mh * m[j]; + // hpk = h[k] - mh * m[k]; + // the micromagnetic standard problem 4 failed to converge (?) + // + // this->NOTE (Fri 08 Jul 2016 13:58): In fact, the problem converges but with 2 less + // decimals of accuracy, compared with the OOMMF calculation + double mth0 = 0, mth1 = 0, mth2 = 0; + + // The first term: m x H_eff = m x H_perp + // if (do_precession){ + if (std::abs(llg_params->gamma) > 0){ // Correct with EPS value + mth0 = cross_x(m[i], m[j], m[k], hpi, hpj, hpk); + mth1 = cross_y(m[i], m[j], m[k], hpi, hpj, hpk); + mth2 = cross_z(m[i], m[j], m[k], hpi, hpj, hpk); + } + + // The RHS term of the LLG equation + dmdt[i] = coeff * (mth0 - hpi * llg_params->alpha[id]); + dmdt[j] = coeff * (mth1 - hpj * llg_params->alpha[id]); + dmdt[k] = coeff * (mth2 - hpk * llg_params->alpha[id]); + + // In future, we will try the new method to integrate the LLG equation, + // A mixed mid-point Runge-Kutta like scheme for the integration of + // Landau-Lifshitz equation Journal of Applied Physics 115, 17D101 + // (2014) if possible, we can combine it with adaptive step size, don't + // know how to do but it's worth a try. + + // if (default_c < 0){ + // c = 6 * sqrt(dm_dt[i] * dm_dt[i] + + // dm_dt[j] * dm_dt[j] + + // dm_dt[k] * dm_dt[k] + // ); + // } else { + // c = default_c; + // } + //printf("%0.15g %0.15g\n", c, default_c); + + // TODO: Correct the RHS term to keep m normalised + c = 1.; + dmdt[i] += c * (1 - mm) * m[i]; + dmdt[j] += c * (1 - mm) * m[j]; + dmdt[k] += c * (1 - mm) * m[k]; + + } + +} + +// Passes the required parameters for the integrator to the LLG_params struct +void MicroLLGDriver::setup(MicroSim * sim, + double * alpha, double gamma, double t, double dt) { + + this->llg_params = new LLG_params; // initialize struct + this->sim = sim; + this->llg_params->gamma = gamma; + this->llg_params->alpha = alpha; + this->llg_params->field = sim->field; + this->llg_params->pins = sim->pins; + this->t = t; // Should be handled by integrator + this->dt = dt; // should be handled by integrator + this->dmdt.resize(3 * this->sim->n); + + // DEBUG + std::cout << "[driver] Gamma:" << this->llg_params->gamma << std::endl; + std::cout << "[driver] Sim n:" << this->sim->n << std::endl; +} + + +void MicroLLGDriver::add_integrator(IntegratorID integrator_id) { + + switch(integrator_id) { + case RK4: + // How to call Derived class integrator if the pointer is to the Base class? + this->integrator = new IntegratorRK4(this->sim->n); + // Best create a set_integrator_parameters function + this->integrator->setup(this->sim->n, this->dt); + // TODO: separate function to set time? + this->integrator->t = 0.0; + this->integrator->step_n = 0; + + std::cout << "Setup integrator" << std::endl; + + } + +} + +void MicroLLGDriver::run_until(double t_final){ + + while(this->integrator->t <= t_final) { + this->integrator->integration_step(LLG_RHS, + this->sim->spin, + this->dmdt, + this->llg_params); + } + +} + +// ---------------------------------------------------------------------------- + +template +void IntegratorRK4::setup(unsigned int N, double dt) { + // Initial values + this->step_n = 0; + this->t = 0.0; + this->dt = dt; + this->N = N; + // this->compute_RHS = f; + + // Initial values?? + // for (unsigned int i = 0; i < this->N; ++i ) { + // rksteps[i] = init_m[i]; + // } +} + +template +void IntegratorRK4::integration_step(void (*f) (double * m, + std::vector& dmdt, + unsigned int n, + double t, + T * params), + double * m, + std::vector& dmdt, + T * params) { + + // To be defined in constructor/setup: + double t_factor[4] = {0.0, 0.5, 0.5, 1.0}; + double m_factor[4] = {0.0, 0.5, 0.5, 1.0}; + std::vector& rksteps = this->integratorData; + + for (unsigned int RKSTEP = 0; RKSTEP < 4; ++RKSTEP) { + + // Re-use the rk_steps vector to apply the RK step + // For RKSTEP = 0 we initialise the first this->N elements of rk_steps vector with m[i] + // Should we separate the step calculation for every i value?? + for (unsigned int i = 0; i < this->N; ++i ) { + unsigned int RKSTEP_idx = i + (RKSTEP) * this->N; + + rksteps[RKSTEP_idx] = m[i]; + + if (RKSTEP > 0) { + rksteps[RKSTEP_idx] += m_factor[RKSTEP] * this->dt * rksteps[i + (RKSTEP - 1) * this->N]; + } + + } + // Update the corresponding values of the RK step vector + f(&rksteps[RKSTEP * this->N], dmdt, this->N, this->t + t_factor[RKSTEP] * this->dt, params); + // Copy RHS values to the corresponding RK step array section + for (unsigned int i = 0; i < this->N; ++i ) rksteps[i + RKSTEP * this->N] = dmdt[i]; + } + + // Update time and array + this->t += this->dt; + for (unsigned int i = 0; i < this->N; ++i) { + m[i] += (1 / 6) * this->dt * ( rksteps[i] + + 2 * rksteps[i + this->N] + + 2 * rksteps[i + 2 * this->N] + + rksteps[i + 3 * this->N]); + } + + // // RK4 step 1 + // this->compute_RHS(rk_steps.begin(), rk_steps.begin() + this->N, rk_steps, t); + + // // RK4 step 2 + // for (unsigned int i = 0; i < this->N; ++i ) { + // rk_steps[i + this->N] = m[i] + 0.5 * this->dt * rk_steps[i]; + // } + // this->compute_RHS(rk_steps.begin() + this->N, rk_steps.begin() + 2 * this->N, + // rk_steps, t + 0.5 * dt); + + // // RK4 step 3 + // for (unsigned int i = 0; i < this->N; ++i ) { + // rk_steps[i + 2 * this->N] = m[i] + 0.5 * this->dt * rk_steps[this->N + i]; + // } + // this->compute_RHS(rk_steps.begin() + 2 * this->N, rk_steps.begin() + 3 * this->N, + // rk_steps, t + 0.5 * dt); + + // // RK4 step 4 + // for (unsigned int i = 0; i < this->N; ++i ) { + // rk_steps[i + 2 * this->N] = m[i] + 0.5 * this->dt * rk_steps[2 * this->N + i]; + // } + // this->compute_RHS(rk_steps.begin() + 3 * this->N, rk_steps.begin() + 4 * this->N, + // rk_steps, t + 1.0 * dt); + +} + +// For testing: --------------------------------------------------------------- + +template class IntegratorRK4; + +// TODO: If we move the integrator to a separate file, we should declare all +// the templates used by Fidimag here, e.g. +// template class IntegratorRK4; diff --git a/native/src/m_driver_test.cpp b/native/src/m_driver_test.cpp new file mode 100644 index 00000000..337114b6 --- /dev/null +++ b/native/src/m_driver_test.cpp @@ -0,0 +1,39 @@ +#include "m_driver_test.h" +#include "c_vectormath.h" +#include + +void test_f(double * x, std::vector& dxdt, unsigned int n, + double t, TestParams * params) { + + // Function: x' = x ((2 / (e^t + 1)) - 1) + // Solution: x(t) = 12 * e^t / (e^t + 1)^2 + for (int i = 0; i < n; ++i) { + dxdt[i] = x[i] * ((2 / (std::exp(t) + 1)) - 1); + } + +} + +void test_integrator_RK4 (void) { + + int N = 10; + double dt = 1e-3; + TestParams * testp; + testp->foo = 1.; + // TestParams * testp = &(TestParams) {.foo = 1.}; -> ? + double x[N]; + std::vector dxdt(N); + + // This particular template is declared in m_driver.cpp for now + IntegratorRK4 integratorrk4(N); + integratorrk4.setup(N, dt); + + // For now set initial m manually: + for (int i = 0; i < N; ++i) { + x[i] = 1.0; + } + + int int_steps = 20; + for (int i = 0; i < N; ++i) { + integratorrk4.integration_step(test_f, x, dxdt, testp); + } +} diff --git a/fidimag/micro/lib/exch.c b/native/src/m_exch.cpp similarity index 54% rename from fidimag/micro/lib/exch.c rename to native/src/m_exch.cpp index a12c4d69..1f6564e8 100644 --- a/fidimag/micro/lib/exch.c +++ b/native/src/m_exch.cpp @@ -1,8 +1,14 @@ -#include "micro_clib.h" +#include "m_clib.h" +// New energy class: +#include "c_energy.h" +#include "c_constants.h" -void compute_exch_field_micro(double *restrict m, double *restrict field, double *restrict energy, - double *restrict Ms_inv, double A, double dx, double dy, double dz, - int n, int *restrict ngbs) { +// ---------------------------------------------------------------------------- +// WILL BE DEPRECATED SOON :: + +void compute_exch_field_micro(double * m, double * field, double * energy, + double * Ms_inv, double A, double dx, double dy, double dz, + int n, int * ngbs) { /* Compute the micromagnetic exchange field and energy using the * matrix of neighbouring spins and a second order approximation @@ -231,6 +237,166 @@ void compute_exch_field_rkky_micro(double *m, double *field, double *energy, dou } } } +} +// ---------------------------------------------------------------------------- +// Using the new energy class:: +void ExchangeEnergy::compute_field(double t) { + /* Compute the micromagnetic exchange field and energy using the + * matrix of neighbouring spins and a second order approximation + * for the derivative + * + * Ms_inv :: Array with the (1 / Ms) values for every mesh node. + * The values are zero for points with Ms = 0 (no material) + * + * A :: Exchange constant + * + * dx, dy, dz :: Mesh spacings in the corresponding directions + * + * n :: Number of mesh nodes + * + * ngbs :: The array of neighbouring spins, which has (6 * n) + * entries. Specifically, it contains the indexes of + * the neighbours of every mesh node, in the following order: + * -x, +x, -y, +y, -z, +z + * + * Thus, the array is like: + * | 0-x, 0+x, 0-y, 0+y, 0-z, 0+z, 1-x, 1+x, 1-y, ... | + * i=0 i=1 ... + * + * where 0-y is the index of the neighbour of the 0th spin, + * in the -y direction, for example. The index value for a + * neighbour where Ms = 0, is evaluated as -1. The array + * automatically gives periodic boundaries. + * + * A basic example is a 3 x 3 two dimensional mesh with PBCs + * in the X and Y direction: + * + * +-----------+ + * | 6 | 7 | 8 | + * +-----------+ + * | 3 | 4 | 5 | + * +-----------+ + * | 0 | 1 | 2 | + * +-----------+ + * + * so, the first 6 entries (neighbours of the 0th mesh node) + * of the array would be: [ 2 1 6 3 -1 -1 ... ] + * (-1 since there is no material in +-z, and a '2' first, + * since it is the left neighbour which is the PBC in x, etc..) + * + * For the exchange computation, the field is defined as: + * H_ex = (2 * A / (mu0 * Ms)) * nabla^2 (mx, my, mz) + * + * Therefore, for the i-th mesh node (spin), we approximate the + * derivatives as: + * nabla^2 mx = (1 / dx^2) * ( spin[i-x] - 2 * spin[i] + spin[i+x] ) + + * (1 / dy^2) * ( spin[i-y] - 2 * spin[i] + spin[i+y] ) + + * (1 / dz^2) * ( spin[i-z] - 2 * spin[i] + spin[i+z] ) + * + * Where i-x is the neighbour in the -x direction. This is similar + * for my and mz. + * We can notice that the sum is the same if we do: + * ( spin[i-x] - spin[i] ) + ( spin[i+x] - spin[i] ) + * so we can iterate through the neighbours and perform the sum with the + * corresponding coefficient 1 /dx, 1/dy or 1/dz + * + * The *m array contains the spins as: + * [mx0, my0, mz0, mx1, my1, mz1, mx2, ...] + * so if we want the starting position of the magnetisation for the + * i-th spin, we only have to do (3 * i) for mx, (3 * i + 1) for my + * and (3 * i + 2) for mz + * + * + * IMPORTANT: The ex field usually has the structure: + * 2 * A / (mu0 Ms ) * (Second derivative of M) + * When discretising the derivative, it carries a "2" in the + * denominator which we "cancel" with the "2" in the prefactor, + * hence we do not put it explicitly in the calculations + * + * So, when computing the energy: (-1/2) * mu * Ms * H_ex + * we only put the 0.5 factor and don't worry about the "2"s in the + * field + * + */ + + /* Here we iterate through every mesh node */ + for (int i = 0; i < n; i++) { + /* Define the coefficients */ + double ax = 2 * A[i] / (dx * dx); + double ay = 2 * A[i] / (dy * dy); + double az = 2 * A[i] / (dz * dz); + + double fx = 0, fy = 0, fz = 0; + int idnm = 0; // Index for the magnetisation matrix + int idn = 6 * i; // index for the neighbours + + /* Set a zero field for sites without magnetic material */ + if (Ms_inv[i] == 0.0){ + field[3 * i] = 0; + field[3 * i + 1] = 0; + field[3 * i + 2] = 0; + continue; + } + + /* Here we iterate through the neighbours */ + for (int j = 0; j < 6; j++) { + /* Remember that index=-1 is for sites without material */ + if (ngbs[idn + j] >= 0) { + /* Magnetisation of the neighbouring spin since ngbs gives + * the neighbour's index */ + idnm = 3 * ngbs[idn + j]; + + /* Check that the magnetisation of the neighbouring spin + * is larger than zero */ + if (Ms_inv[ngbs[idn + j]] > 0){ + + /* Neighbours in the -x and +x directions + * giving: ( spin[i-x] - spin[i] ) + ( spin[i+x] - spin[i] ) + * when ngbs[idn + j] > 0 for j = 0 and j=1 + * If, for example, there is no + * neighbour at -x (j=0) in the 0th node (no PBCs), + * the second derivative would only be avaluated as: + * (1 / dx * dx) * ( spin[i+x] - spin[i] ) + * which, according to + * [M.J. Donahue and D.G. Porter; Physica B, 343, 177-183 (2004)] + * when performing the integration of the energy, we still + * have error of the order O(dx^2) + * This same applies for the other directions + */ + if (j == 0 || j == 1) { + fx += ax * (spin[idnm] - spin[3 * i]); + fy += ax * (spin[idnm + 1] - spin[3 * i + 1]); + fz += ax * (spin[idnm + 2] - spin[3 * i + 2]); + } + /* Neighbours in the -y and +y directions */ + else if (j == 2 || j == 3) { + fx += ay * (spin[idnm] - spin[3 * i]); + fy += ay * (spin[idnm + 1] - spin[3 * i + 1]); + fz += ay * (spin[idnm + 2] - spin[3 * i + 2]); + } + /* Neighbours in the -z and +z directions */ + else if (j == 4 || j == 5) { + fx += az * (spin[idnm] - spin[3 * i]); + fy += az * (spin[idnm + 1] - spin[3 * i + 1]); + fz += az * (spin[idnm + 2] - spin[3 * i + 2]); + } + else { + continue; + } + } + } + } + + /* Energy as: (-mu0 * Ms / 2) * [ H_ex * m ] */ + energy[i] = -0.5 * (fx * spin[3 * i] + fy * spin[3 * i + 1] + + fz * spin[3 * i + 2]); + + /* Update the field H_ex which has the same structure than *m */ + field[3 * i] = fx * Ms_inv[i] * MU_0_INV; + field[3 * i + 1] = fy * Ms_inv[i] * MU_0_INV; + field[3 * i + 2] = fz * Ms_inv[i] * MU_0_INV; + + } } diff --git a/fidimag/micro/lib/util.c b/native/src/m_util.cpp similarity index 90% rename from fidimag/micro/lib/util.c rename to native/src/m_util.cpp index 6982fe3e..362d4116 100644 --- a/fidimag/micro/lib/util.c +++ b/native/src/m_util.cpp @@ -1,4 +1,6 @@ -#include "micro_clib.h" +#include "m_clib.h" +#include "c_vectormath.h" +#include "c_constants.h" // Compute: S \cdot (S_i \times S_j) inline double volume(double S[3], double Si[3], double Sj[3]) { @@ -8,8 +10,8 @@ inline double volume(double S[3], double Si[3], double Sj[3]) { return tx + ty + tz; } -double skyrmion_number(double *restrict spin, double *restrict charge, - int nx, int ny, int nz, int *restrict ngbs) { +double skyrmion_number(double * spin, double * charge, + int nx, int ny, int nz, int * ngbs) { /* Compute the skyrmion number Q, defined as: * _ @@ -35,7 +37,7 @@ double skyrmion_number(double *restrict spin, double *restrict charge, * * The *spin array is the vector field for a two dimensional * lattice with dimensions nx * ny - * (we can take a slice of a bulk from Python and pass it here, + * (we can take a slice of a bulk from Python and pass it here, * remember to do the ame for the neighbours matrix) * The array follows the order: * [Sx0 Sy0 Sz0 Sx1 Sy1 Sz1 ... ] @@ -63,9 +65,9 @@ double skyrmion_number(double *restrict spin, double *restrict charge, /* Store the spin directions of the nearest neighbours * in the order: [-x +x -y +y] - */ + */ double S_nn[12]; - + for(i=0;i<12;i++){ S_nn[i] = 0; //we have to set S_nn to zeros manually } @@ -73,7 +75,7 @@ double skyrmion_number(double *restrict spin, double *restrict charge, for (i = 0; i < nxy; i++) { index = 3 * i; - /* The starting index of the nearest neighbours for the + /* The starting index of the nearest neighbours for the * i-th spin */ int id_nn = 6 * i; @@ -95,7 +97,7 @@ double skyrmion_number(double *restrict spin, double *restrict charge, - volume(S, &S_nn[0], &S_nn[9]) - volume(S, &S_nn[3], &S_nn[6]); - charge[i] /= (16 * WIDE_PI); + charge[i] /= (16 * M_PIl); sum += charge[i]; } diff --git a/setup.py b/setup.py index 4bf2ef69..df09ba49 100644 --- a/setup.py +++ b/setup.py @@ -1,259 +1,60 @@ from distutils.core import setup from distutils.extension import Extension import multiprocessing -# from Cython.Distutils import build_ext +from tools import glob_files, BuildError, get_user_module_sources from Cython.Build import cythonize import numpy -import os import glob -import re -import sys -#if sys.platform == 'darwin': -# from distutils import sysconfig -# vars = sysconfig.get_config_vars() -# vars['LDSHARED'] = vars['LDSHARED'].replace('-bundle', '-dynamiclib') - -class BuildError(Exception): - pass +import os -if 'CC' in os.environ: - print("Using CC={}".format(os.environ['CC'])) -else: - os.environ["CC"] = "gcc" - print("Using CC={} (set by setup.py)".format(os.environ['CC'])) +print(os.environ['CC'], os.environ['CXX']) +version = '5.0alpha' MODULE_DIR = os.path.dirname(os.path.abspath(__file__)) -print(MODULE_DIR) -SRC_DIR = os.path.join(MODULE_DIR, "fidimag") -SUNDIALS_DIR = os.path.join(SRC_DIR, "common", "sundials") -NEB_DIR = os.path.join(SRC_DIR, "common", "neb") -NEBM_DIR = os.path.join(SRC_DIR, "common", "neb_method") -ATOM_DIR = os.path.join(SRC_DIR, "atomistic", "lib") -COMMON_DIR = os.path.join(SRC_DIR, "common", "lib") -MICRO_DIR = os.path.join(SRC_DIR, "micro", "lib") -BARYAKHTAR_DIR = os.path.join(MICRO_DIR, "baryakhtar") -DEMAG_DIR = os.path.join(SRC_DIR, "common", "dipolar") -USER_DIR = os.path.join(SRC_DIR, "user") -print(USER_DIR) - -LOCAL_DIR = os.path.join(MODULE_DIR, "local") -INCLUDE_DIR = os.path.join(LOCAL_DIR, "include") -LIB_DIR = os.path.join(LOCAL_DIR, "lib") -print("LIB_DIR={}".format(LIB_DIR)) - - -pkg_init_path = os.path.join( - os.path.dirname(__file__), 'fidimag', '__init__.py') - - -def get_version(): - with open(pkg_init_path) as f: - for line in f: - m = re.match(r'''__version__\s*=\s*(['"])(.+)\1''', line.strip()) - if m: - return m.group(2) - raise Exception("Couldn't find __version__ in %s" % pkg_init_path) - -version = get_version() - - -def glob_cfiles(path, excludes, extension="*.c"): - cfiles = [] - for cfile in glob.glob(os.path.join(path, extension)): - filename = os.path.basename(cfile) - print(filename) - if not filename in tuple(excludes): - cfiles.append(cfile) - return cfiles - -sources = [] -sources.append(os.path.join(ATOM_DIR, 'clib.pyx')) -sources += glob_cfiles(ATOM_DIR, excludes=["clib.c"]) - - - -common_sources = [] -common_sources.append(os.path.join(COMMON_DIR, 'common_clib.pyx')) -common_sources += glob_cfiles(COMMON_DIR, excludes=["common_clib.c"]) - -cvode_sources = [] -cvode_sources.append(os.path.join(SUNDIALS_DIR, 'cvode.pyx')) -cvode_sources += glob_cfiles(SUNDIALS_DIR, excludes=["cvode.c"]) - -baryakhtar_sources = [] -baryakhtar_sources.append(os.path.join(BARYAKHTAR_DIR, 'baryakhtar_clib.pyx')) -baryakhtar_sources += glob_cfiles(BARYAKHTAR_DIR, - excludes=["baryakhtar_clib.c"]) +INCLUDE_DIR = os.path.join(MODULE_DIR, 'local', 'include') +LIB_DIR = os.path.join(MODULE_DIR, 'local', 'lib') -micro_sources = [] -micro_sources.append(os.path.join(MICRO_DIR, 'micro_clib.pyx')) -micro_sources += glob_cfiles(MICRO_DIR, excludes=["micro_clib.c"]) - -# NEB Method ------------------------------------------------------------------ - -nebm_sources = [] -nebm_sources.append(os.path.join(NEBM_DIR, "nebm_clib.pyx")) -nebm_sources += glob_cfiles(NEBM_DIR, excludes=["nebm_clib.c"]) - -# ----------------------------------------------------------------------------- - -dipolar_sources = [] -dipolar_sources.append(os.path.join(DEMAG_DIR, 'dipolar.pyx')) -dipolar_sources += glob_cfiles(DEMAG_DIR, excludes=["dipolar.c"]) - - -com_libs = ['m', 'fftw3_omp', 'fftw3', 'sundials_cvodes', - 'sundials_nvecserial', 'sundials_nvecopenmp', 'blas', 'lapack'] - -com_args = ['-std=c99', '-O3', '-Wno-cpp', '-Wno-unused-function'] -# rpath is the path relative to the compiled shared object files (e.g. clib.so, etc) -# which the dynamic linker looks for the linked libraries (e.g. libsundials_*.so) in. -# We need to set it relatively in order for it to be preserved if the parent directory is moved -# hence why it is a 'relative'(r) path. Here the relative path is with respect to -# the fidimag/fidimag/extensions directory. RPATH = '../../local/lib' com_link = ['-Wl,-rpath,{}'.format(LIB_DIR)] -lib_paths = [LIB_DIR] - -if 'icc' in os.environ['CC']: - com_args.append('-openmp') - com_link.append('-openmp') -else: - com_args.append('-fopenmp') +lib_paths = [LIB_DIR, os.path.join(MODULE_DIR, 'native')] +com_inc = [numpy.get_include(), INCLUDE_DIR, os.path.join(MODULE_DIR, 'native', 'include')] +com_libs = ['fidimag'] +com_args = [] - - - com_link.append('-fopenmp') - - -com_inc = [numpy.get_include(), INCLUDE_DIR] - -if 'SUNDIALS_DIR' in os.environ: - lib_paths.append(os.environ['SUNDIALS_DIR']) - com_inc.append(os.environ['SUNDIALS_INC']) +if 'SUNDIALS_INC' in os.environ: + com_inc.append(os.environ['SUNDIALS_INC']) if 'FFTW_DIR' in os.environ: - lib_paths.append(os.environ['FFTW_DIR']) - com_inc.append(os.environ['FFTW_INC']) - -ext_modules = [ - Extension("fidimag.extensions.clib", - sources=sources, - include_dirs=com_inc, - libraries=com_libs, - library_dirs=lib_paths, runtime_library_dirs=lib_paths, - extra_compile_args=com_args, - extra_link_args=com_link, - ), - Extension("fidimag.extensions.common_clib", - sources=common_sources, - include_dirs=com_inc, - libraries=com_libs, - library_dirs=lib_paths, runtime_library_dirs=lib_paths, - extra_compile_args=com_args, - extra_link_args=com_link, - ), - Extension("fidimag.extensions.cvode", - sources=cvode_sources, - include_dirs=com_inc, - libraries=com_libs, - library_dirs=lib_paths, runtime_library_dirs=lib_paths, - extra_compile_args=com_args, - extra_link_args=com_link, - ), - Extension("fidimag.extensions.baryakhtar_clib", - sources=baryakhtar_sources, - include_dirs=com_inc, - libraries=com_libs, - library_dirs=lib_paths, runtime_library_dirs=lib_paths, - extra_compile_args=com_args, - extra_link_args=com_link, - ), - Extension("fidimag.extensions.micro_clib", - sources=micro_sources, - include_dirs=com_inc, - libraries=com_libs, - library_dirs=lib_paths, runtime_library_dirs=lib_paths, - extra_compile_args=com_args, - extra_link_args=com_link, - ), - Extension("fidimag.extensions.nebm_clib", - sources=nebm_sources, - include_dirs=com_inc, - libraries=com_libs, - library_dirs=lib_paths, runtime_library_dirs=lib_paths, - extra_compile_args=com_args, - extra_link_args=com_link, - ), - Extension("fidimag.extensions.cvode", - sources=cvode_sources, - include_dirs=com_inc, - libraries=com_libs, - library_dirs=lib_paths, runtime_library_dirs=lib_paths, - extra_compile_args=com_args, - extra_link_args=com_link, - ), - Extension("fidimag.extensions.baryakhtar_clib", - sources=baryakhtar_sources, - include_dirs=com_inc, - libraries=com_libs, - library_dirs=lib_paths, runtime_library_dirs=lib_paths, - extra_compile_args=com_args, - extra_link_args=com_link, - ), - Extension("fidimag.extensions.micro_clib", - sources=micro_sources, - include_dirs=com_inc, - libraries=com_libs, - library_dirs=lib_paths, runtime_library_dirs=lib_paths, - extra_compile_args=com_args, - extra_link_args=com_link, - ), - Extension("fidimag.extensions.dipolar", - sources=dipolar_sources, - include_dirs=com_inc, - libraries=com_libs, - library_dirs=lib_paths, runtime_library_dirs=lib_paths, - extra_compile_args=com_args, - extra_link_args=com_link, - ), -] + com_inc.append(os.environ['FFTW_INC']) -for folder in glob.glob(os.path.join(USER_DIR, '*/')): - module_name = folder.split('/')[-2] - print('Found User Module: {}'.format(module_name)) - user_sources = glob.glob(folder + '/*.pyx') - print('\tFound Cython sources: {}'.format(user_sources)) - if len(user_sources) != 1: - raise BuildError("User Modules are only allowed one Cython .pyx file") +source_files = glob.glob(os.path.join('fidimag', '**', '*.pyx'), recursive=True) +print(source_files) +ext_names = ["fidimag.extensions." + s.split('/')[-1].rstrip('.pyx') for s in source_files] - filename_string = user_sources[0].split('/')[-1][:-4] - if filename_string != module_name: - print(filename_string, module_name) - raise BuildError("The Cython source file in {} must match the folder name - i.e. it must be {}.pyx".format(module_name, module_name)) - cfilename = filename_string + '.c' - print(cfilename) - user_sources += glob_cfiles(folder, excludes=[cfilename]) +print(ext_names) - print(user_sources) - ext_modules.append( - Extension("fidimag.extensions.user.{}".format(module_name), - sources=user_sources, - include_dirs=com_inc, - libraries=com_libs, - library_dirs=lib_paths, runtime_library_dirs=lib_paths, - extra_compile_args=com_args, - extra_link_args=com_link, - ), +ext_modules = [] +for i, (module, src) in enumerate(zip(ext_names, source_files)): + print("Compiling module {}".format(module)) + ext_modules.append(Extension(module, + sources=[src], + include_dirs=com_inc, + libraries=com_libs, + library_dirs=lib_paths, runtime_library_dirs=lib_paths, + extra_compile_args=com_args, + extra_link_args=com_link, ) +) + + nthreads = multiprocessing.cpu_count() -print('Building with {} threads'.format(nthreads)) + setup( name='fidimag', version=version, @@ -263,6 +64,12 @@ def glob_cfiles(path, excludes, extension="*.c"): 'fidimag.micro', 'fidimag.extensions', 'fidimag.common', - ], - ext_modules=cythonize(ext_modules, nthreads=nthreads), + ], + ext_modules=cythonize(ext_modules, + nthreads=nthreads, + compiler_directives={ + 'linetrace': True, + 'language_level': '3', + } + ), ) diff --git a/tests/test_C_classes.py b/tests/test_C_classes.py new file mode 100644 index 00000000..08719a26 --- /dev/null +++ b/tests/test_C_classes.py @@ -0,0 +1,55 @@ +# See: fidimag/common/c_clib.pyx for Python wrappers +from fidimag.extensions.c_clib import PyExchangeEnergy +from fidimag.extensions.c_clib import PyMicroSim +from fidimag.extensions.c_clib import PyMicroLLGDriver +from fidimag.extensions.c_clib import INTEGRATOR_RK4 +from fidimag.common import CuboidMesh + +import time +import numpy as np + +nx, ny, nz = 3, 3, 1 +n = nx * ny * nz +dx, dy, dz = 1, 1, 1 +unit_length = 1. + +mesh = CuboidMesh(nx=nx, ny=ny, nz=nz, dx=dx, dy=dy, dz=dz) + +spin = np.ones(3 * n, dtype=np.double) +spin[::3] = 0 +spin[1::3] = 0 +spin[3 * 4:3 * 4 + 3] = [0, 1, 0] +Ms = np.ones(n) +Ms_inv = 1 / Ms + +field = np.zeros(3 * n, dtype=np.double) +energy = np.zeros(n, dtype=np.double) +pins = np.zeros(n, dtype=np.int32) + +sim_C = PyMicroSim() + +sim_C.setup(mesh.nx, mesh.ny, mesh.nz, mesh.dx, mesh.dy, mesh.dz, + mesh.unit_length, mesh.coordinates, mesh.neighbours, + spin, Ms, Ms_inv, energy, field, pins) + +A = np.ones(9, dtype=np.double) +print(A) +# A[2] = 4 +# Exch.printA() + +Exch = PyExchangeEnergy(A, sim_C) +# Exch.setup(sim_C) + +Exch.compute_field(0) +print(field) +print(energy) + +sim_C.add_interaction(Exch) + +# time.sleep(5) + +alpha = np.ones(n) +gamma = 1. +driver = PyMicroLLGDriver() +driver.setup(sim_C, alpha, gamma, t=0.0, dt=0.1) +driver.add_integrator(INTEGRATOR_RK4) diff --git a/tests/test_monte_carlo.py b/tests/test_monte_carlo.py index 5f4c6f58..2f1c1f53 100644 --- a/tests/test_monte_carlo.py +++ b/tests/test_monte_carlo.py @@ -2,7 +2,7 @@ mpl.use("Agg") import matplotlib.pyplot as plt -import fidimag.extensions.clib as clib +import fidimag.extensions.a_clib as clib import numpy as np from fidimag.atomistic import MonteCarlo, HexagonalMesh import fidimag.common.constant as const diff --git a/tests/test_skyrmion_number.py b/tests/test_skyrmion_number.py index 1e9aebc1..611bbb14 100644 --- a/tests/test_skyrmion_number.py +++ b/tests/test_skyrmion_number.py @@ -49,7 +49,7 @@ def init_m_multiple_sks(pos, r, sk_pos): def test_skx_num_atomistic(): - """ + r""" Test the *finite spin chirality* or skyrmion number for a discrete spins simulation in a two dimensional lattice diff --git a/tests/test_stt.py b/tests/test_stt.py index ffa9f593..88ecde98 100644 --- a/tests/test_stt.py +++ b/tests/test_stt.py @@ -1,9 +1,9 @@ import numpy as np -import fidimag.extensions.common_clib as clib +import fidimag.extensions.c_clib as clib def test_sst_field_1d(): - """ + r""" This is a direct test of the STT C library We create a 1-D 4 spins system along the x direction with the following components: diff --git a/tools.py b/tools.py new file mode 100644 index 00000000..0c18e7bc --- /dev/null +++ b/tools.py @@ -0,0 +1,35 @@ +import glob +import os + +def glob_files(path, excludes, extension="*.cpp"): + files = [] + for srcfile in glob.glob(os.path.join(path, extension)): + filename = os.path.basename(srcfile) + if not filename in tuple(excludes): + files.append(srcfile) + return files + + +class BuildError(Exception): + pass + + + +def get_user_module_sources(folder): + module_name = folder.split('/')[-2] + print('Found User Module: {}'.format(folder)) + user_sources = glob.glob(folder + '/*.pyx') + print('\tFound Cython sources: {}'.format(user_sources)) + + if len(user_sources) != 1: + raise BuildError("User Modules are only allowed one Cython .pyx file") + + filename_string = user_sources[0].split('/')[-1][:-4] + if filename_string != module_name: + print(filename_string, module_name) + raise BuildError("The Cython source file in {} must match the folder name - i.e. it must be {}.pyx".format(module_name, module_name)) + cfilename = filename_string + '.cpp' + print(cfilename) + user_sources += glob_files(folder, excludes=[cfilename]) + + return module_name, user_sources