diff --git a/fidimag/common/dipolar/demag.c b/fidimag/common/dipolar/demag.c index 1c908ae2..0dfd5707 100644 --- a/fidimag/common/dipolar/demag.c +++ b/fidimag/common/dipolar/demag.c @@ -1,512 +1,485 @@ +#include "demagcoef.h" +#include "dipolar.h" #include -#include #include -#include "dipolar.h" -#include "demagcoef.h" - +#include inline double Nxxdipole(double x, double y, double z) { - double x2 = x * x; - double y2 = y * y; - double z2 = z * z; - double R = x2 + y2 + z2; - if (R == 0) - return 0.0; - double r = sqrt(R); - return -(2 * x2 - y2 - z2) / (R * R * r); + double x2 = x * x; + double y2 = y * y; + double z2 = z * z; + double R = x2 + y2 + z2; + if (R == 0) + return 0.0; + double r = sqrt(R); + return -(2 * x2 - y2 - z2) / (R * R * r); } inline double Nxydipole(double x, double y, double z) { - double R = x * x + y * y + z * z; - if (R == 0) - return 0.0; - double r = sqrt(R); - return -3 * x * y / (R * R * r); + double R = x * x + y * y + z * z; + if (R == 0) + return 0.0; + double r = sqrt(R); + return -3 * x * y / (R * R * r); } double NXXdipole(enum Type_Nij type, double x, double y, double z) { - switch (type) { - case Tensor_xx: - return Nxxdipole(x, y, z); - case Tensor_yy: - return Nxxdipole(y, x, z); - case Tensor_zz: - return Nxxdipole(z, y, x); - case Tensor_xy: - return Nxydipole(x, y, z); - case Tensor_xz: - return Nxydipole(x, z, y); - case Tensor_yz: - return Nxydipole(y, z, x); - } - return 0; + switch (type) { + case Tensor_xx: + return Nxxdipole(x, y, z); + case Tensor_yy: + return Nxxdipole(y, x, z); + case Tensor_zz: + return Nxxdipole(z, y, x); + case Tensor_xy: + return Nxydipole(x, y, z); + case Tensor_xz: + return Nxydipole(x, z, y); + case Tensor_yz: + return Nxydipole(y, z, x); + } + return 0; } -//compute the demag tensors, i.e, H=-N.M +// compute the demag tensors, i.e, H=-N.M void compute_dipolar_tensors(fft_demag_plan *plan) { - int i, j, k, id; - double x, y, z; - - int nx = plan->nx; - int ny = plan->ny; - int nz = plan->nz; - int lenx = plan->lenx; - int leny = plan->leny; - int lenz = plan->lenz; - int lenxy = lenx * leny; - double dx = plan->dx; - double dy = plan->dy; - double dz = plan->dz; - // Parallelising this like this - // means that z should be the largest index - // in order to get better performance from threading. - // The data writing is not very clever here; we're - // going to be invalidating the cache a lot. It would be better - // maybe to split this into six seperate loops for each component - // of the tensor in order that each thread is working on a smaller - // memory address range - #pragma omp parallel for private(j, i, x, y, z, id) schedule(dynamic, 32) - for (k = 0; k < lenz; k++) { - for (j = 0; j < leny; j++) { - for (i = 0; i < lenx; i++) { - id = k * lenxy + j * lenx + i; - - if ((lenx%2 == 0 && i == nx) || (leny%2 == 0 && j == ny) || (lenz%2 == 0 && k == nz)) { - plan->tensor_xx[id] = 0.0; - plan->tensor_yy[id] = 0.0; - plan->tensor_zz[id] = 0.0; - plan->tensor_xy[id] = 0.0; - plan->tensor_xz[id] = 0.0; - plan->tensor_yz[id] = 0.0; - continue; - } - - x = (itensor_xx[id] = NXXdipole(Tensor_xx, x, y, z); - plan->tensor_yy[id] = NXXdipole(Tensor_yy, x, y, z); - plan->tensor_zz[id] = NXXdipole(Tensor_zz, x, y, z); - plan->tensor_xy[id] = NXXdipole(Tensor_xy, x, y, z); - plan->tensor_xz[id] = NXXdipole(Tensor_xz, x, y, z); - plan->tensor_yz[id] = NXXdipole(Tensor_yz, x, y, z); - - } - } - } + int i, j, k, id; + double x, y, z; + + int nx = plan->nx; + int ny = plan->ny; + int nz = plan->nz; + int lenx = plan->lenx; + int leny = plan->leny; + int lenz = plan->lenz; + int lenxy = lenx * leny; + double dx = plan->dx; + double dy = plan->dy; + double dz = plan->dz; +// Parallelising this like this +// means that z should be the largest index +// in order to get better performance from threading. +// The data writing is not very clever here; we're +// going to be invalidating the cache a lot. It would be better +// maybe to split this into six seperate loops for each component +// of the tensor in order that each thread is working on a smaller +// memory address range +#pragma omp parallel for private(j, i, x, y, z, id) schedule(dynamic, 32) + for (k = 0; k < lenz; k++) { + for (j = 0; j < leny; j++) { + for (i = 0; i < lenx; i++) { + id = k * lenxy + j * lenx + i; + + if ((lenx % 2 == 0 && i == nx) || (leny % 2 == 0 && j == ny) || (lenz % 2 == 0 && k == nz)) { + plan->tensor_xx[id] = 0.0; + plan->tensor_yy[id] = 0.0; + plan->tensor_zz[id] = 0.0; + plan->tensor_xy[id] = 0.0; + plan->tensor_xz[id] = 0.0; + plan->tensor_yz[id] = 0.0; + continue; + } + + x = (i < nx) ? i * dx : (i - lenx) * dx; + y = (j < ny) ? j * dy : (j - leny) * dy; + z = (k < nz) ? k * dz : (k - lenz) * dz; + + plan->tensor_xx[id] = NXXdipole(Tensor_xx, x, y, z); + plan->tensor_yy[id] = NXXdipole(Tensor_yy, x, y, z); + plan->tensor_zz[id] = NXXdipole(Tensor_zz, x, y, z); + plan->tensor_xy[id] = NXXdipole(Tensor_xy, x, y, z); + plan->tensor_xz[id] = NXXdipole(Tensor_xz, x, y, z); + plan->tensor_yz[id] = NXXdipole(Tensor_yz, x, y, z); + } + } + } } -//compute the demag tensors, i.e, H=-N.M +// compute the demag tensors, i.e, H=-N.M void compute_demag_tensors(fft_demag_plan *plan) { - int i, j, k, id; - double x, y, z, radius_sq; - - int nx = plan->nx; - int ny = plan->ny; - int nz = plan->nz; - int lenx = plan->lenx; - int leny = plan->leny; - int lenz = plan->lenz; - int lenxy = lenx * leny; - - double dx = plan->dx; - double dy = plan->dy; - double dz = plan->dz; - - double length = pow(dx*dy*dz, 1/3.0); - double asymptotic_radius_sq = pow(26.0*length,2.0); - #pragma omp parallel for private(j, i, id, x, y, z, radius_sq) schedule(dynamic, 32) - for (k = 0; k < lenz; k++) { - for (j = 0; j < leny; j++) { - for (i = 0; i < lenx; i++) { - id = k * lenxy + j * lenx + i; - - if ((lenx%2 == 0 && i == nx) || (leny%2 == 0 && j == ny) || (lenz%2 == 0 && k == nz)) { - plan->tensor_xx[id] = 0.0; - plan->tensor_yy[id] = 0.0; - plan->tensor_zz[id] = 0.0; - plan->tensor_xy[id] = 0.0; - plan->tensor_xz[id] = 0.0; - plan->tensor_yz[id] = 0.0; - continue; - } - - x = (iasymptotic_radius_sq){ - //printf("%g %g %g %g %g %g\n",x,y,z,dx,dy,dz); - plan->tensor_xx[id] = DemagNxxAsymptotic(x, y, z, dx, dy, dz); - plan->tensor_yy[id] = DemagNyyAsymptotic(x, y, z, dx, dy, dz); - plan->tensor_zz[id] = DemagNzzAsymptotic(x, y, z, dx, dy, dz); - plan->tensor_xy[id] = DemagNxyAsymptotic(x, y, z, dx, dy, dz); - plan->tensor_xz[id] = DemagNxzAsymptotic(x, y, z, dx, dy, dz); - plan->tensor_yz[id] = DemagNyzAsymptotic(x, y, z, dx, dy, dz); - }else{ - //printf("%g %g %g %g %g %g\n",x,y,z,dx,dy,dz); - plan->tensor_xx[id] = CalculateSDA00(x, y, z, dx, dy, dz); - plan->tensor_yy[id] = CalculateSDA11(x, y, z, dx, dy, dz); - plan->tensor_zz[id] = CalculateSDA22(x, y, z, dx, dy, dz); - plan->tensor_xy[id] = CalculateSDA01(x, y, z, dx, dy, dz); - plan->tensor_xz[id] = CalculateSDA02(x, y, z, dx, dy, dz); - plan->tensor_yz[id] = CalculateSDA12(x, y, z, dx, dy, dz); - } - } - } - } + int i, j, k, id; + double x, y, z, radius_sq; + + int nx = plan->nx; + int ny = plan->ny; + int nz = plan->nz; + int lenx = plan->lenx; + int leny = plan->leny; + int lenz = plan->lenz; + int lenxy = lenx * leny; + + double dx = plan->dx; + double dy = plan->dy; + double dz = plan->dz; + + double length = pow(dx * dy * dz, 1 / 3.0); + double asymptotic_radius_sq = pow(26.0 * length, 2.0); +#pragma omp parallel for private(j, i, id, x, y, z, radius_sq) schedule(dynamic, 32) + for (k = 0; k < lenz; k++) { + for (j = 0; j < leny; j++) { + for (i = 0; i < lenx; i++) { + id = k * lenxy + j * lenx + i; + + if ((lenx % 2 == 0 && i == nx) || (leny % 2 == 0 && j == ny) || (lenz % 2 == 0 && k == nz)) { + plan->tensor_xx[id] = 0.0; + plan->tensor_yy[id] = 0.0; + plan->tensor_zz[id] = 0.0; + plan->tensor_xy[id] = 0.0; + plan->tensor_xz[id] = 0.0; + plan->tensor_yz[id] = 0.0; + continue; + } + + x = (i < nx) ? i * dx : (i - lenx) * dx; + y = (j < ny) ? j * dy : (j - leny) * dy; + z = (k < nz) ? k * dz : (k - lenz) * dz; + + radius_sq = x * x + y * y + z * z; + + if (radius_sq > asymptotic_radius_sq) { + // printf("%g %g %g %g %g %g\n",x,y,z,dx,dy,dz); + plan->tensor_xx[id] = DemagNxxAsymptotic(x, y, z, dx, dy, dz); + plan->tensor_yy[id] = DemagNyyAsymptotic(x, y, z, dx, dy, dz); + plan->tensor_zz[id] = DemagNzzAsymptotic(x, y, z, dx, dy, dz); + plan->tensor_xy[id] = DemagNxyAsymptotic(x, y, z, dx, dy, dz); + plan->tensor_xz[id] = DemagNxzAsymptotic(x, y, z, dx, dy, dz); + plan->tensor_yz[id] = DemagNyzAsymptotic(x, y, z, dx, dy, dz); + } else { + // printf("%g %g %g %g %g %g\n",x,y,z,dx,dy,dz); + plan->tensor_xx[id] = CalculateSDA00(x, y, z, dx, dy, dz); + plan->tensor_yy[id] = CalculateSDA11(x, y, z, dx, dy, dz); + plan->tensor_zz[id] = CalculateSDA22(x, y, z, dx, dy, dz); + plan->tensor_xy[id] = CalculateSDA01(x, y, z, dx, dy, dz); + plan->tensor_xz[id] = CalculateSDA02(x, y, z, dx, dy, dz); + plan->tensor_yz[id] = CalculateSDA12(x, y, z, dx, dy, dz); + } + } + } + } } - - -//used for debug +// used for debug void print_r(char *str, double *restrict x, int n) { - int i; - printf("%s:\n", str); - for (i = 0; i < n; i++) { - printf("%f ", x[i]); - } - printf("\n"); - + int i; + printf("%s:\n", str); + for (i = 0; i < n; i++) { + printf("%f ", x[i]); + } + printf("\n"); } void print_c(char *str, fftw_complex *restrict x, int n) { - int i; - printf("%s\n", str); - for (i = 0; i < n; i++) { - //printf("%f+%fI ", x[i]); - } - printf("\n"); - + int i; + printf("%s\n", str); + for (i = 0; i < n; i++) { + // printf("%f+%fI ", x[i]); + } + printf("\n"); } fft_demag_plan *create_plan(void) { - fft_demag_plan *plan = (fft_demag_plan*) malloc(sizeof(fft_demag_plan)); + fft_demag_plan *plan = (fft_demag_plan *)malloc(sizeof(fft_demag_plan)); - return plan; + return plan; } void init_plan(fft_demag_plan *restrict plan, double dx, double dy, - double dz, int nx, int ny, int nz) { + double dz, int nx, int ny, int nz) { - //plan->mu_s = mu_s; + // plan->mu_s = mu_s; - fftw_init_threads(); - fftw_plan_with_nthreads(omp_get_max_threads()); + fftw_init_threads(); + fftw_plan_with_nthreads(omp_get_max_threads()); - plan->dx = dx; - plan->dy = dy; - plan->dz = dz; + plan->dx = dx; + plan->dy = dy; + plan->dz = dz; - plan->nx = nx; - plan->ny = ny; - plan->nz = nz; - - int critical_n = 8; + plan->nx = nx; + plan->ny = ny; + plan->nz = nz; - plan->lenx = nx > critical_n ? 2 * nx : 2 * nx - 1; - plan->leny = ny > critical_n ? 2 * ny : 2 * ny - 1; - plan->lenz = nz > critical_n ? 2 * nz : 2 * nz - 1; + int critical_n = 8; - plan->total_length = plan->lenx * plan->leny * plan->lenz; + plan->lenx = nx > critical_n ? 2 * nx : 2 * nx - 1; + plan->leny = ny > critical_n ? 2 * ny : 2 * ny - 1; + plan->lenz = nz > critical_n ? 2 * nz : 2 * nz - 1; - int size1 = plan->total_length * sizeof(double); - int size2 = plan->total_length * sizeof(fftw_complex); + plan->total_length = plan->lenx * plan->leny * plan->lenz; - plan->tensor_xx = (double *) fftw_malloc(size1); - plan->tensor_yy = (double *) fftw_malloc(size1); - plan->tensor_zz = (double *) fftw_malloc(size1); - plan->tensor_xy = (double *) fftw_malloc(size1); - plan->tensor_xz = (double *) fftw_malloc(size1); - plan->tensor_yz = (double *) fftw_malloc(size1); + int size1 = plan->total_length * sizeof(double); + int size2 = plan->total_length * sizeof(fftw_complex); - plan->mx = (double *) fftw_malloc(size1); - plan->my = (double *) fftw_malloc(size1); - plan->mz = (double *) fftw_malloc(size1); + plan->tensor_xx = (double *)fftw_malloc(size1); + plan->tensor_yy = (double *)fftw_malloc(size1); + plan->tensor_zz = (double *)fftw_malloc(size1); + plan->tensor_xy = (double *)fftw_malloc(size1); + plan->tensor_xz = (double *)fftw_malloc(size1); + plan->tensor_yz = (double *)fftw_malloc(size1); - plan->hx = (double *) fftw_malloc(size1); - plan->hy = (double *) fftw_malloc(size1); - plan->hz = (double *) fftw_malloc(size1); + plan->mx = (double *)fftw_malloc(size1); + plan->my = (double *)fftw_malloc(size1); + plan->mz = (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->hx = (double *)fftw_malloc(size1); + plan->hy = (double *)fftw_malloc(size1); + plan->hz = (double *)fftw_malloc(size1); - 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->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->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); } - void create_fftw_plan(fft_demag_plan *restrict 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->m_plan = fftw_plan_dft_r2c_3d(plan->lenz, plan->leny, plan->lenx, - plan->mx, 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); - - for (int i = 0; i < plan->total_length; i++) { - plan->Nxx[i] = 0; - plan->Nyy[i] = 0; - plan->Nzz[i] = 0; - plan->Nxy[i] = 0; - plan->Nxz[i] = 0; - plan->Nyz[i] = 0; - - plan->mx[i] = 0; - plan->my[i] = 0; - plan->mz[i] = 0; - - plan->hx[i] = 0; - plan->hy[i] = 0; - plan->hz[i] = 0; - } - - - 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_destroy_plan(plan->tensor_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->m_plan = fftw_plan_dft_r2c_3d(plan->lenz, plan->leny, plan->lenx, + plan->mx, 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); + + for (int i = 0; i < plan->total_length; i++) { + plan->Nxx[i] = 0; + plan->Nyy[i] = 0; + plan->Nzz[i] = 0; + plan->Nxy[i] = 0; + plan->Nxz[i] = 0; + plan->Nyz[i] = 0; + + plan->mx[i] = 0; + plan->my[i] = 0; + plan->mz[i] = 0; + + plan->hx[i] = 0; + plan->hy[i] = 0; + plan->hz[i] = 0; + } + + 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_destroy_plan(plan->tensor_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 +// 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) { - int i, j, k, id1, id2; - - int nx = plan->nx; - int ny = plan->ny; - int nz = plan->nz; - int nxy = nx * ny; - - int lenx = plan->lenx; - int leny = plan->leny; - int lenxy = lenx * leny; - - for (i = 0; i < plan->total_length; i++) { - plan->mx[i] = 0; - plan->my[i] = 0; - plan->mz[i] = 0; - } - - - for (k = 0; k < nz; k++) { - for (j = 0; j < ny; j++) { - for (i = 0; i < nx; i++) { - id1 = k * nxy + j * nx + i; - id2 = k * lenxy + j * lenx + i; - - plan->mx[id2] = spin[3*id1]*mu_s[id1]; - plan->my[id2] = spin[3*id1+1]*mu_s[id1]; - plan->mz[id2] = spin[3*id1+2]*mu_s[id1]; - } - } - } - - //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); - - //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; - - //print_c("Mx", Mx, plan->total_length); - - for (i = 0; i < plan->total_length; i++) { - Hx[i] = Nxx[i] * Mx[i] + Nxy[i] * My[i] + Nxz[i] * Mz[i]; - Hy[i] = Nxy[i] * Mx[i] + Nyy[i] * My[i] + Nyz[i] * Mz[i]; - Hz[i] = Nxz[i] * Mx[i] + Nyz[i] * My[i] + Nzz[i] * Mz[i]; - } - - //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); - //print_r("hx", plan->hx, plan->total_length); - //print_r("hy", plan->hy, plan->total_length); - //print_r("hz", plan->hz, plan->total_length); - - double scale = -1.0 / plan->total_length; - #pragma omp parallel for private(j, i, id1, id2) schedule(dynamic, 32) - for (k = 0; k < nz; k++) { - for (j = 0; j < ny; j++) { - for (i = 0; i < nx; i++) { - id1 = k * nxy + j * nx + i; - - id2 = k * lenxy + j * lenx + i; - field[3*id1] = plan->hx[id2] * scale; - field[3*id1+1] = plan->hy[id2] * scale; - field[3*id1+2] = plan->hz[id2] * scale; - } - } - } - + int i, j, k, id1, id2; + + int nx = plan->nx; + int ny = plan->ny; + int nz = plan->nz; + int nxy = nx * ny; + + int lenx = plan->lenx; + int leny = plan->leny; + int lenxy = lenx * leny; + + for (i = 0; i < plan->total_length; i++) { + plan->mx[i] = 0; + plan->my[i] = 0; + plan->mz[i] = 0; + } + + for (k = 0; k < nz; k++) { + for (j = 0; j < ny; j++) { + for (i = 0; i < nx; i++) { + id1 = k * nxy + j * nx + i; + id2 = k * lenxy + j * lenx + i; + + plan->mx[id2] = spin[3 * id1] * mu_s[id1]; + plan->my[id2] = spin[3 * id1 + 1] * mu_s[id1]; + plan->mz[id2] = spin[3 * id1 + 2] * mu_s[id1]; + } + } + } + + // 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); + + // 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; + + // print_c("Mx", Mx, plan->total_length); + + for (i = 0; i < plan->total_length; i++) { + Hx[i] = Nxx[i] * Mx[i] + Nxy[i] * My[i] + Nxz[i] * Mz[i]; + Hy[i] = Nxy[i] * Mx[i] + Nyy[i] * My[i] + Nyz[i] * Mz[i]; + Hz[i] = Nxz[i] * Mx[i] + Nyz[i] * My[i] + Nzz[i] * Mz[i]; + } + + // 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); + // print_r("hx", plan->hx, plan->total_length); + // print_r("hy", plan->hy, plan->total_length); + // print_r("hz", plan->hz, plan->total_length); + + double scale = -1.0 / plan->total_length; +#pragma omp parallel for private(j, i, id1, id2) schedule(dynamic, 32) + for (k = 0; k < nz; k++) { + for (j = 0; j < ny; j++) { + for (i = 0; i < nx; i++) { + id1 = k * nxy + j * nx + i; + + id2 = k * lenxy + j * lenx + i; + field[3 * id1] = plan->hx[id2] * scale; + field[3 * id1 + 1] = plan->hy[id2] * scale; + field[3 * id1 + 2] = plan->hz[id2] * scale; + } + } + } } -//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; - - - 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); - } - } - } - +// 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; + + 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); + } + } + } } double compute_demag_energy(fft_demag_plan *restrict plan, double *restrict spin, double *restrict mu_s, double *restrict field, - double *restrict energy - ) { - - int i,j; + double *restrict energy) { - int nxyz = plan->nx * plan->ny * plan->nz; + int i, j; - double total_energy = 0; + int nxyz = plan->nx * plan->ny * plan->nz; - for (i = 0; i < nxyz; i++) { - j = 3*i; - energy[i] = -0.5 * mu_s[i] * (spin[j] * field[j] + - spin[j+1] * field[j+1] + - spin[j+2] * field[j+2] ); - total_energy += energy[i]; - } + double total_energy = 0; - return total_energy; + for (i = 0; i < nxyz; i++) { + j = 3 * i; + energy[i] = -0.5 * mu_s[i] * (spin[j] * field[j] + spin[j + 1] * field[j + 1] + spin[j + 2] * field[j + 2]); + total_energy += energy[i]; + } + return total_energy; } void finalize_plan(fft_demag_plan *restrict plan) { - fftw_destroy_plan(plan->m_plan); - fftw_destroy_plan(plan->h_plan); - - fftw_free(plan->tensor_xx); - fftw_free(plan->tensor_yy); - fftw_free(plan->tensor_zz); - fftw_free(plan->tensor_xy); - fftw_free(plan->tensor_xz); - fftw_free(plan->tensor_yz); - - fftw_free(plan->Nxx); - fftw_free(plan->Nyy); - fftw_free(plan->Nzz); - fftw_free(plan->Nxy); - fftw_free(plan->Nxz); - fftw_free(plan->Nyz); - - fftw_free(plan->Mx); - fftw_free(plan->My); - fftw_free(plan->Mz); - fftw_free(plan->Hx); - fftw_free(plan->Hy); - fftw_free(plan->Hz); - - fftw_free(plan->mx); - fftw_free(plan->my); - fftw_free(plan->mz); - fftw_free(plan->hx); - fftw_free(plan->hy); - fftw_free(plan->hz); - - fftw_cleanup_threads(); - - free(plan); + fftw_destroy_plan(plan->m_plan); + fftw_destroy_plan(plan->h_plan); + + fftw_free(plan->tensor_xx); + fftw_free(plan->tensor_yy); + fftw_free(plan->tensor_zz); + fftw_free(plan->tensor_xy); + fftw_free(plan->tensor_xz); + fftw_free(plan->tensor_yz); + + fftw_free(plan->Nxx); + fftw_free(plan->Nyy); + fftw_free(plan->Nzz); + fftw_free(plan->Nxy); + fftw_free(plan->Nxz); + fftw_free(plan->Nyz); + + fftw_free(plan->Mx); + fftw_free(plan->My); + fftw_free(plan->Mz); + fftw_free(plan->Hx); + fftw_free(plan->Hy); + fftw_free(plan->Hz); + + fftw_free(plan->mx); + fftw_free(plan->my); + fftw_free(plan->mz); + fftw_free(plan->hx); + fftw_free(plan->hy); + fftw_free(plan->hz); + + fftw_cleanup_threads(); + + free(plan); } - diff --git a/fidimag/common/dipolar/demag_oommf.c b/fidimag/common/dipolar/demag_oommf.c index 8675802d..deb25824 100644 --- a/fidimag/common/dipolar/demag_oommf.c +++ b/fidimag/common/dipolar/demag_oommf.c @@ -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 @@ -39,12 +39,10 @@ to copyright protection within the United States. */ - +#include "demagcoef.h" +#include "dipolar.h" #include #include -#include "dipolar.h" -#include "demagcoef.h" - /* FILE: demagcoef.cc -*-Mode: c++-*- * @@ -120,119 +118,117 @@ to copyright protection within the United States. * */ - //////////////////////////////////////////////////////////////////////////// // Routines to calculate kernel coefficients // See Newell et al. for details. The code below follows the // naming conventions in that paper. -double SelfDemagNx(double x,double y,double z) -{ // Here Hx = -Nxx.Mx (formula (16) in Newell). +double SelfDemagNx(double x, double y, double z) { // Here Hx = -Nxx.Mx (formula (16) in Newell). // Note: egcs-2.91.57 on Linux/x86 with -O1 mangles this // function (produces NaN's) unless we manually group terms. - if(x<=0.0 || y<=0.0 || z<=0.0) + if (x <= 0.0 || y <= 0.0 || z <= 0.0) return 0.0; - if (x==y && y==z) - return 1./3.; // Special case: cube + if (x == y && y == z) + return 1. / 3.; // Special case: cube - double xsq=x*x,ysq=y*y,zsq=z*z; - double diag=sqrt(xsq+ysq+zsq); + double xsq = x * x, ysq = y * y, zsq = z * z; + double diag = sqrt(xsq + ysq + zsq); double arr[15]; - double mpxy = (x-y)*(x+y); - double mpxz = (x-z)*(x+z); + double mpxy = (x - y) * (x + y); + double mpxz = (x - z) * (x + z); - arr[0] = -4*(2*xsq*x-ysq*y-zsq*z); - arr[1] = 4*(xsq+mpxy)*sqrt(xsq+ysq); - arr[2] = 4*(xsq+mpxz)*sqrt(xsq+zsq); - arr[3] = -4*(ysq+zsq)*sqrt(ysq+zsq); - arr[4] = -4*diag*(mpxy+mpxz); + arr[0] = -4 * (2 * xsq * x - ysq * y - zsq * z); + arr[1] = 4 * (xsq + mpxy) * sqrt(xsq + ysq); + arr[2] = 4 * (xsq + mpxz) * sqrt(xsq + zsq); + arr[3] = -4 * (ysq + zsq) * sqrt(ysq + zsq); + arr[4] = -4 * diag * (mpxy + mpxz); - arr[5] = 24*x*y*z*atan(y*z/(x*diag)); - arr[6] = 12*(z+y)*xsq*log(x); + arr[5] = 24 * x * y * z * atan(y * z / (x * diag)); + arr[6] = 12 * (z + y) * xsq * log(x); - arr[7] = 12*z*ysq*log((sqrt(ysq+zsq)+z)/y); - arr[8] = -12*z*xsq*log(sqrt(xsq+zsq)+z); - arr[9] = 12*z*mpxy*log(diag+z); - arr[10] = -6*z*mpxy*log(xsq+ysq); + arr[7] = 12 * z * ysq * log((sqrt(ysq + zsq) + z) / y); + arr[8] = -12 * z * xsq * log(sqrt(xsq + zsq) + z); + arr[9] = 12 * z * mpxy * log(diag + z); + arr[10] = -6 * z * mpxy * log(xsq + ysq); - arr[11] = 12*y*zsq*log((sqrt(ysq+zsq)+y)/z); - arr[12] = -12*y*xsq*log(sqrt(xsq+ysq)+y); - arr[13] = 12*y*mpxz*log(diag+y); - arr[14] = -6*y*mpxz*log(xsq+zsq); + arr[11] = 12 * y * zsq * log((sqrt(ysq + zsq) + y) / z); + arr[12] = -12 * y * xsq * log(sqrt(xsq + ysq) + y); + arr[13] = 12 * y * mpxz * log(diag + y); + arr[14] = -6 * y * mpxz * log(xsq + zsq); - double Nxx = AccurateSum(15,arr)/(12*WIDE_PI*x*y*z); + double Nxx = AccurateSum(15, arr) / (12 * WIDE_PI * x * y * z); return Nxx; } -double SelfDemagNy(double xsize,double ysize,double zsize) -{ return SelfDemagNx(ysize,zsize,xsize); } - -double SelfDemagNz(double xsize,double ysize,double zsize) -{ return SelfDemagNx(zsize,xsize,ysize); } +double SelfDemagNy(double xsize, double ysize, double zsize) { return SelfDemagNx(ysize, zsize, xsize); } +double SelfDemagNz(double xsize, double ysize, double zsize) { return SelfDemagNx(zsize, xsize, ysize); } double -Newell_f(double x,double y,double z) -{ // There is mucking around here to handle case where imports +Newell_f(double x, double y, double z) { // There is mucking around here to handle case where imports // are near zero. In particular, asinh(t) is written as // log(t+sqrt(1+t)) because the latter appears easier to // handle if t=y/x (for example) as x -> 0. - // This function is even; the fabs()'s just simplify special case handling. - x=fabs(x); double xsq=x*x; - y=fabs(y); double ysq=y*y; - z=fabs(z); double zsq=z*z; + // This function is even; the fabs()'s just simplify special case handling. + x = fabs(x); + double xsq = x * x; + y = fabs(y); + double ysq = y * y; + z = fabs(z); + double zsq = z * z; - double R=xsq+ysq+zsq; - if(R<=0.0) return 0.0; - else R=sqrt(R); + double R = xsq + ysq + zsq; + if (R <= 0.0) + return 0.0; + else + R = sqrt(R); // f(x,y,z) double piece[8]; - int piececount=0; - if(z>0.) { // For 2D grids, half the calls from F1 have z==0. - double temp1,temp2,temp3; - piece[piececount++] = 2*(2*xsq-ysq-zsq)*R; - if((temp1=x*y*z)>0.) - piece[piececount++] = -12*temp1*atan2(y*z,x*R); - if(y>0. && (temp2=xsq+zsq)>0.) { - double dummy = log(((y+R)*(y+R))/temp2); - piece[piececount++] = 3*y*zsq*dummy; - piece[piececount++] = -3*y*xsq*dummy; + int piececount = 0; + if (z > 0.) { // For 2D grids, half the calls from F1 have z==0. + double temp1, temp2, temp3; + piece[piececount++] = 2 * (2 * xsq - ysq - zsq) * R; + if ((temp1 = x * y * z) > 0.) + piece[piececount++] = -12 * temp1 * atan2(y * z, x * R); + if (y > 0. && (temp2 = xsq + zsq) > 0.) { + double dummy = log(((y + R) * (y + R)) / temp2); + piece[piececount++] = 3 * y * zsq * dummy; + piece[piececount++] = -3 * y * xsq * dummy; } - if((temp3=xsq+ysq)>0.) { - double dummy = log(((z+R)*(z+R))/temp3); - piece[piececount++] = 3*z*ysq*dummy; - piece[piececount++] = -3*z*xsq*dummy; + if ((temp3 = xsq + ysq) > 0.) { + double dummy = log(((z + R) * (z + R)) / temp3); + piece[piececount++] = 3 * z * ysq * dummy; + piece[piececount++] = -3 * z * xsq * dummy; } } else { // z==0 - if(x==y) { - //const double K = 2*sqrt(static_cast(2.0)) - // -6*log(1+sqrt(static_cast(2.0))); + if (x == y) { + // const double K = 2*sqrt(static_cast(2.0)) + // -6*log(1+sqrt(static_cast(2.0))); double K = -2.4598143973710680537922785014593576970294; - piece[piececount++] = K*xsq*x; + piece[piececount++] = K * xsq * x; } else { - piece[piececount++] = 2*(2*xsq-ysq)*R; - if(y>0. && x>0.) - piece[piececount++] = -6*y*xsq*log((y+R)/x); + piece[piececount++] = 2 * (2 * xsq - ysq) * R; + if (y > 0. && x > 0.) + piece[piececount++] = -6 * y * xsq * log((y + R) / x); } } - return AccurateSum(piececount,piece)/12.; + return AccurateSum(piececount, piece) / 12.; } double -CalculateSDA00(double x,double y,double z, - double dx,double dy,double dz) -{ // This is Nxx in Newell's paper - double result=0.; - if(x==0. && y==0. && z==0.) { +CalculateSDA00(double x, double y, double z, + double dx, double dy, double dz) { // This is Nxx in Newell's paper + double result = 0.; + if (x == 0. && y == 0. && z == 0.) { // Self demag term. The base routine can handle x==y==z==0, // but this should be more accurate. - result = SelfDemagNx(dx,dy,dz); + result = SelfDemagNx(dx, dy, dz); } else { // Simplified (collapsed) formula based on Newell's paper. // This saves about half the calls to f(). There is still @@ -240,96 +236,98 @@ CalculateSDA00(double x,double y,double z, // but as this is an initialization-only issue speed shouldn't // be too critical. double arr[27]; - arr[ 0] = -1*Newell_f(x+dx,y+dy,z+dz); - arr[ 1] = -1*Newell_f(x+dx,y-dy,z+dz); - arr[ 2] = -1*Newell_f(x+dx,y-dy,z-dz); - arr[ 3] = -1*Newell_f(x+dx,y+dy,z-dz); - arr[ 4] = -1*Newell_f(x-dx,y+dy,z-dz); - arr[ 5] = -1*Newell_f(x-dx,y+dy,z+dz); - arr[ 6] = -1*Newell_f(x-dx,y-dy,z+dz); - arr[ 7] = -1*Newell_f(x-dx,y-dy,z-dz); - - arr[ 8] = 2*Newell_f(x,y-dy,z-dz); - arr[ 9] = 2*Newell_f(x,y-dy,z+dz); - arr[10] = 2*Newell_f(x,y+dy,z+dz); - arr[11] = 2*Newell_f(x,y+dy,z-dz); - arr[12] = 2*Newell_f(x+dx,y+dy,z); - arr[13] = 2*Newell_f(x+dx,y,z+dz); - arr[14] = 2*Newell_f(x+dx,y,z-dz); - arr[15] = 2*Newell_f(x+dx,y-dy,z); - arr[16] = 2*Newell_f(x-dx,y-dy,z); - arr[17] = 2*Newell_f(x-dx,y,z+dz); - arr[18] = 2*Newell_f(x-dx,y,z-dz); - arr[19] = 2*Newell_f(x-dx,y+dy,z); - - arr[20] = -4*Newell_f(x,y-dy,z); - arr[21] = -4*Newell_f(x,y+dy,z); - arr[22] = -4*Newell_f(x,y,z-dz); - arr[23] = -4*Newell_f(x,y,z+dz); - arr[24] = -4*Newell_f(x+dx,y,z); - arr[25] = -4*Newell_f(x-dx,y,z); - - arr[26] = 8*Newell_f(x,y,z); - - result = AccurateSum(27,arr)/(4*WIDE_PI*dx*dy*dz); + arr[0] = -1 * Newell_f(x + dx, y + dy, z + dz); + arr[1] = -1 * Newell_f(x + dx, y - dy, z + dz); + arr[2] = -1 * Newell_f(x + dx, y - dy, z - dz); + arr[3] = -1 * Newell_f(x + dx, y + dy, z - dz); + arr[4] = -1 * Newell_f(x - dx, y + dy, z - dz); + arr[5] = -1 * Newell_f(x - dx, y + dy, z + dz); + arr[6] = -1 * Newell_f(x - dx, y - dy, z + dz); + arr[7] = -1 * Newell_f(x - dx, y - dy, z - dz); + + arr[8] = 2 * Newell_f(x, y - dy, z - dz); + arr[9] = 2 * Newell_f(x, y - dy, z + dz); + arr[10] = 2 * Newell_f(x, y + dy, z + dz); + arr[11] = 2 * Newell_f(x, y + dy, z - dz); + arr[12] = 2 * Newell_f(x + dx, y + dy, z); + arr[13] = 2 * Newell_f(x + dx, y, z + dz); + arr[14] = 2 * Newell_f(x + dx, y, z - dz); + arr[15] = 2 * Newell_f(x + dx, y - dy, z); + arr[16] = 2 * Newell_f(x - dx, y - dy, z); + arr[17] = 2 * Newell_f(x - dx, y, z + dz); + arr[18] = 2 * Newell_f(x - dx, y, z - dz); + arr[19] = 2 * Newell_f(x - dx, y + dy, z); + + arr[20] = -4 * Newell_f(x, y - dy, z); + arr[21] = -4 * Newell_f(x, y + dy, z); + arr[22] = -4 * Newell_f(x, y, z - dz); + arr[23] = -4 * Newell_f(x, y, z + dz); + arr[24] = -4 * Newell_f(x + dx, y, z); + arr[25] = -4 * Newell_f(x - dx, y, z); + + arr[26] = 8 * Newell_f(x, y, z); + + result = AccurateSum(27, arr) / (4 * WIDE_PI * dx * dy * dz); } return result; } - double -Newell_g(double x,double y,double z) -{ // There is mucking around here to handle case where imports +Newell_g(double x, double y, double z) { // There is mucking around here to handle case where imports // are near zero. In particular, asinh(t) is written as // log(t+sqrt(1+t)) because the latter appears easier to // handle if t=y/x (for example) as x -> 0. - double result_sign=1.0; - if(x<0.0) - result_sign *= -1.0; - if(y<0.0) + double result_sign = 1.0; + if (x < 0.0) result_sign *= -1.0; - x=fabs(x); y=fabs(y); z=fabs(z); // This function is even in z and + 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) + double xsq = x * x, ysq = y * y, zsq = z * z; + double R = xsq + ysq + zsq; + if (R <= 0.0) return 0.0; else - R=sqrt(R); + R = sqrt(R); // g(x,y,z) double piece[7]; - int piececount=0; - piece[piececount++] = -2*x*y*R;; - if(z>0.) { // For 2D grids, 1/3 of the calls from CalculateSDA01 have z==0. - piece[piececount++] = -z*zsq*atan2(x*y,z*R); - piece[piececount++] = -3*z*ysq*atan2(x*z,y*R); - piece[piececount++] = -3*z*xsq*atan2(y*z,x*R); + int piececount = 0; + piece[piececount++] = -2 * x * y * R; + ; + if (z > 0.) { // For 2D grids, 1/3 of the calls from CalculateSDA01 have z==0. + piece[piececount++] = -z * zsq * atan2(x * y, z * R); + piece[piececount++] = -3 * z * ysq * atan2(x * z, y * R); + piece[piececount++] = -3 * z * xsq * atan2(y * z, x * R); - double temp1,temp2,temp3; - if((temp1=xsq+ysq)>0.) - piece[piececount++] = 6*x*y*z*log((z+R)/sqrt(temp1)); + double temp1, temp2, temp3; + if ((temp1 = xsq + ysq) > 0.) + piece[piececount++] = 6 * x * y * z * log((z + R) / sqrt(temp1)); - if((temp2=ysq+zsq)>0.) - piece[piececount++] = y*(3*zsq-ysq)*log((x+R)/sqrt(temp2)); + if ((temp2 = ysq + zsq) > 0.) + piece[piececount++] = y * (3 * zsq - ysq) * log((x + R) / sqrt(temp2)); - if((temp3=xsq+zsq)>0.) - piece[piececount++] = x*(3*zsq-xsq)*log((y+R)/sqrt(temp3)); + if ((temp3 = xsq + zsq) > 0.) + piece[piececount++] = x * (3 * zsq - xsq) * log((y + R) / sqrt(temp3)); } else { // z==0. - if(y>0.) piece[piececount++] = -y*ysq*log((x+R)/y); - if(x>0.) piece[piececount++] = -x*xsq*log((y+R)/x); + if (y > 0.) + piece[piececount++] = -y * ysq * log((x + R) / y); + if (x > 0.) + piece[piececount++] = -x * xsq * log((y + R) / x); } - return result_sign*AccurateSum(piececount,piece)/6.; + return result_sign * AccurateSum(piececount, piece) / 6.; } -double CalculateSDA01(double x,double y,double z, - double l,double h,double e) -{ // This is Nxy*(4*PI*tau) in Newell's paper. +double CalculateSDA01(double x, double y, double z, + double l, double h, double e) { // This is Nxy*(4*PI*tau) in Newell's paper. // Simplified (collapsed) formula based on Newell's paper. // This saves about half the calls to g(). There is still @@ -338,200 +336,190 @@ double CalculateSDA01(double x,double y,double z, // be too critical. double arr[27]; - arr[ 0] = -1*Newell_g(x-l,y-h,z-e); - arr[ 1] = -1*Newell_g(x-l,y-h,z+e); - arr[ 2] = -1*Newell_g(x+l,y-h,z+e); - arr[ 3] = -1*Newell_g(x+l,y-h,z-e); - arr[ 4] = -1*Newell_g(x+l,y+h,z-e); - arr[ 5] = -1*Newell_g(x+l,y+h,z+e); - arr[ 6] = -1*Newell_g(x-l,y+h,z+e); - arr[ 7] = -1*Newell_g(x-l,y+h,z-e); - - arr[ 8] = 2*Newell_g(x,y+h,z-e); - arr[ 9] = 2*Newell_g(x,y+h,z+e); - arr[10] = 2*Newell_g(x,y-h,z+e); - arr[11] = 2*Newell_g(x,y-h,z-e); - arr[12] = 2*Newell_g(x-l,y-h,z); - arr[13] = 2*Newell_g(x-l,y+h,z); - arr[14] = 2*Newell_g(x-l,y,z-e); - arr[15] = 2*Newell_g(x-l,y,z+e); - arr[16] = 2*Newell_g(x+l,y,z+e); - arr[17] = 2*Newell_g(x+l,y,z-e); - arr[18] = 2*Newell_g(x+l,y-h,z); - arr[19] = 2*Newell_g(x+l,y+h,z); - - arr[20] = -4*Newell_g(x-l,y,z); - arr[21] = -4*Newell_g(x+l,y,z); - arr[22] = -4*Newell_g(x,y,z+e); - arr[23] = -4*Newell_g(x,y,z-e); - arr[24] = -4*Newell_g(x,y-h,z); - arr[25] = -4*Newell_g(x,y+h,z); - - arr[26] = 8*Newell_g(x,y,z); - - return AccurateSum(27,arr)/(4*WIDE_PI*l*h*e); + arr[0] = -1 * Newell_g(x - l, y - h, z - e); + arr[1] = -1 * Newell_g(x - l, y - h, z + e); + arr[2] = -1 * Newell_g(x + l, y - h, z + e); + arr[3] = -1 * Newell_g(x + l, y - h, z - e); + arr[4] = -1 * Newell_g(x + l, y + h, z - e); + arr[5] = -1 * Newell_g(x + l, y + h, z + e); + arr[6] = -1 * Newell_g(x - l, y + h, z + e); + arr[7] = -1 * Newell_g(x - l, y + h, z - e); + + arr[8] = 2 * Newell_g(x, y + h, z - e); + arr[9] = 2 * Newell_g(x, y + h, z + e); + arr[10] = 2 * Newell_g(x, y - h, z + e); + arr[11] = 2 * Newell_g(x, y - h, z - e); + arr[12] = 2 * Newell_g(x - l, y - h, z); + arr[13] = 2 * Newell_g(x - l, y + h, z); + arr[14] = 2 * Newell_g(x - l, y, z - e); + arr[15] = 2 * Newell_g(x - l, y, z + e); + arr[16] = 2 * Newell_g(x + l, y, z + e); + arr[17] = 2 * Newell_g(x + l, y, z - e); + arr[18] = 2 * Newell_g(x + l, y - h, z); + arr[19] = 2 * Newell_g(x + l, y + h, z); + + arr[20] = -4 * Newell_g(x - l, y, z); + arr[21] = -4 * Newell_g(x + l, y, z); + arr[22] = -4 * Newell_g(x, y, z + e); + arr[23] = -4 * Newell_g(x, y, z - e); + arr[24] = -4 * Newell_g(x, y - h, z); + arr[25] = -4 * Newell_g(x, y + h, z); + + arr[26] = 8 * Newell_g(x, y, z); + + return AccurateSum(27, arr) / (4 * WIDE_PI * l * h * e); /// factor Nxy from Newell's paper. } double DemagNxxAsymptotic(double x, double y, double z, - double dx,double dy,double dz) -{ // Asymptotic approximation to Nxx term. - double R2 = x*x + y*y + z*z; + double dx, double dy, double dz) { // Asymptotic approximation to Nxx term. + double R2 = x * x + y * y + z * z; - if(R2<=0.0) { + if (R2 <= 0.0) { // Asymptotic expansion doesn't apply for R==0. Fall back // to self-demag calculation. - return SelfDemagNx(dx,dy,dz); + return SelfDemagNx(dx, dy, dz); } - double R = sqrt(R2); + double R = sqrt(R2); - double sx2 = x*x/R2; - double sy2 = y*y/R2; - double sz2 = z*z/R2; + double sx2 = x * x / R2; + double sy2 = y * y / R2; + double sz2 = z * z / R2; - double sx4 = sx2*sx2; - double sy4 = sy2*sy2; - double sz4 = sz2*sz2; + double sx4 = sx2 * sx2; + double sy4 = sy2 * sy2; + double sz4 = sz2 * sz2; - double sx6 = sx4*sx2; - double sy6 = sy4*sy2; - double sz6 = sz4*sz2; + double sx6 = sx4 * sx2; + double sy6 = sy4 * sy2; + double sz6 = sz4 * sz2; - double dx2 = dx*dx; - double dy2 = dy*dy; - double dz2 = dz*dz; + double dx2 = dx * dx; + double dy2 = dy * dy; + double dz2 = dz * dz; - double dx4 = dx2*dx2; - double dy4 = dy2*dy2; - double dz4 = dz2*dz2; - double dx2dy2 = dx2*dy2; - double dy2dz2 = dy2*dz2; - double dx2dz2 = dx2*dz2; + double dx4 = dx2 * dx2; + double dy4 = dy2 * dy2; + double dz4 = dz2 * dz2; + double dx2dy2 = dx2 * dy2; + double dy2dz2 = dy2 * dz2; + double dx2dz2 = dx2 * dz2; - double term3 = 2*sx2 - sy2 - sz2; + double term3 = 2 * sx2 - sy2 - sz2; double term5 = 0.0; - if(dx2!=dy2 || dx2!=dz2 || dy2!=dz2) { // Non-cube case - double a1 = 8*dx2 - 4*dy2 - 4*dz2; - double a2 = -24*dx2 + 27*dy2 - 3*dz2; - double a3 = -24*dx2 - 3*dy2 + 27*dz2; - double a4 = 3*dx2 - 4*dy2 + 1*dz2; - double a5 = 6*dx2 - 3*dy2 - 3*dz2; - double a6 = 3*dx2 + 1*dy2 - 4*dz2; - term5 = a1*sx4 + a2*sx2*sy2 + a3*sx2*sz2 - + a4*sy4 + a5*sy2*sz2 + a6*sz4; + if (dx2 != dy2 || dx2 != dz2 || dy2 != dz2) { // Non-cube case + double a1 = 8 * dx2 - 4 * dy2 - 4 * dz2; + double a2 = -24 * dx2 + 27 * dy2 - 3 * dz2; + double a3 = -24 * dx2 - 3 * dy2 + 27 * dz2; + double a4 = 3 * dx2 - 4 * dy2 + 1 * dz2; + double a5 = 6 * dx2 - 3 * dy2 - 3 * dz2; + double a6 = 3 * dx2 + 1 * dy2 - 4 * dz2; + term5 = a1 * sx4 + a2 * sx2 * sy2 + a3 * sx2 * sz2 + a4 * sy4 + a5 * sy2 * sz2 + a6 * sz4; 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; - double b3 = -240*dx4 + 20*dx2dy2 + 580*dx2dz2 + 22*dy4 - 75*dy2dz2 - 202*dz4; - double b4 = 180*dx4 - 505*dx2dy2 + 55*dx2dz2 + 232*dy4 - 75*dy2dz2 + 8*dz4; - double b5 = 360*dx4 - 450*dx2dy2 - 450*dx2dz2 - 180*dy4 + 900*dy2dz2 - 180*dz4; - double b6 = 180*dx4 + 55*dx2dy2 - 505*dx2dz2 + 8*dy4 - 75*dy2dz2 + 232*dz4; - double b7 = -10*dx4 + 30*dx2dy2 - 5*dx2dz2 - 16*dy4 + 10*dy2dz2 - 2*dz4; - double b8 = -30*dx4 + 55*dx2dy2 + 20*dx2dz2 + 8*dy4 - 75*dy2dz2 + 22*dz4; - double b9 = -30*dx4 + 20*dx2dy2 + 55*dx2dz2 + 22*dy4 - 75*dy2dz2 + 8*dz4; - double b10 = -10*dx4 - 5*dx2dy2 + 30*dx2dz2 - 2*dy4 + 10*dy2dz2 - 16*dz4; - - term7 = b1*sx6 + b2*sx4*sy2 + b3*sx4*sz2 + b4*sx2*sy4 + b5*sx2*sy2*sz2 - + b6*sx2*sz4 + b7*sy6 + b8*sy4*sz2 + b9*sy2*sz4 + b10*sz6; - term7 *= (1.0/16.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; + double b3 = -240 * dx4 + 20 * dx2dy2 + 580 * dx2dz2 + 22 * dy4 - 75 * dy2dz2 - 202 * dz4; + double b4 = 180 * dx4 - 505 * dx2dy2 + 55 * dx2dz2 + 232 * dy4 - 75 * dy2dz2 + 8 * dz4; + double b5 = 360 * dx4 - 450 * dx2dy2 - 450 * dx2dz2 - 180 * dy4 + 900 * dy2dz2 - 180 * dz4; + double b6 = 180 * dx4 + 55 * dx2dy2 - 505 * dx2dz2 + 8 * dy4 - 75 * dy2dz2 + 232 * dz4; + double b7 = -10 * dx4 + 30 * dx2dy2 - 5 * dx2dz2 - 16 * dy4 + 10 * dy2dz2 - 2 * dz4; + double b8 = -30 * dx4 + 55 * dx2dy2 + 20 * dx2dz2 + 8 * dy4 - 75 * dy2dz2 + 22 * dz4; + double b9 = -30 * dx4 + 20 * dx2dy2 + 55 * dx2dz2 + 22 * dy4 - 75 * dy2dz2 + 8 * dz4; + double b10 = -10 * dx4 - 5 * dx2dy2 + 30 * dx2dz2 - 2 * dy4 + 10 * dy2dz2 - 16 * dz4; + + term7 = b1 * sx6 + b2 * sx4 * sy2 + b3 * sx4 * sz2 + b4 * sx2 * sy4 + b5 * sx2 * sy2 * sz2 + b6 * sx2 * sz4 + b7 * sy6 + b8 * sy4 * sz2 + b9 * sy2 * sz4 + b10 * sz6; + 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 * WIDE_PI)) * (((term7 / R2 + term5) / R2 + term3) / (R2 * R)); // Error should be of order 1/R^9 return Nxx; } double DemagNxyAsymptotic(double x, double y, double z, - double dx,double dy,double dz) -{ // Asympotic approximation to Nxy term. - double R2 = x*x + y*y + z*z; + double dx, double dy, double dz) { // Asympotic approximation to Nxy term. + double R2 = x * x + y * y + z * z; - if(R2<=0.0) { + if (R2 <= 0.0) { // Asymptotic expansion doesn't apply for R==0. Fall back // to self-demag calculation. return (0.0); } - double R = sqrt(R2); + double R = sqrt(R2); - double sx2 = x*x/R2; - double sy2 = y*y/R2; - double sz2 = z*z/R2; + double sx2 = x * x / R2; + double sy2 = y * y / R2; + double sz2 = z * z / R2; - double sx4 = sx2*sx2; - double sy4 = sy2*sy2; - double sz4 = sz2*sz2; + double sx4 = sx2 * sx2; + double sy4 = sy2 * sy2; + double sz4 = sz2 * sz2; - double dx2 = dx*dx; - double dy2 = dy*dy; - double dz2 = dz*dz; + double dx2 = dx * dx; + double dy2 = dy * dy; + double dz2 = dz * dz; - double dx4 = dx2*dx2; - double dy4 = dy2*dy2; - double dz4 = dz2*dz2; + double dx4 = dx2 * dx2; + double dy4 = dy2 * dy2; + double dz4 = dz2 * dz2; - double dx2dy2 = dx2*dy2; - double dy2dz2 = dy2*dz2; - double dx2dz2 = dx2*dz2; + double dx2dy2 = dx2 * dy2; + double dy2dz2 = dy2 * dz2; + double dx2dz2 = dx2 * dz2; double term3 = 3; double term5 = 0.0; - if(dx2!=dy2 || dx2!=dz2 || dy2!=dz2) { // Non-cube case - double a1 = 4*dx2 - 3*dy2 - 1*dz2; - double a2 = -3*dx2 + 4*dy2 - 1*dz2; - double a3 = -3*dx2 - 3*dy2 + 6*dz2; - term5 = a1*sx2 + a2*sy2 + a3*sz2; - term5 *= (5.0/4.0); + if (dx2 != dy2 || dx2 != dz2 || dy2 != dz2) { // Non-cube case + double a1 = 4 * dx2 - 3 * dy2 - 1 * dz2; + double a2 = -3 * dx2 + 4 * dy2 - 1 * dz2; + double a3 = -3 * dx2 - 3 * dy2 + 6 * dz2; + term5 = a1 * sx2 + a2 * sy2 + a3 * sz2; + term5 *= (5.0 / 4.0); } double term7 = 0.0; { - double b1 = 16*dx4 - 30*dx2dy2 - 10*dx2dz2 + 10*dy4 + 5*dy2dz2 + 2*dz4; - double b2 = -40*dx4 + 105*dx2dy2 - 5*dx2dz2 - 40*dy4 - 5*dy2dz2 + 4*dz4; - double b3 = -40*dx4 - 15*dx2dy2 + 115*dx2dz2 + 20*dy4 - 35*dy2dz2 - 32*dz4; - double b4 = 10*dx4 - 30*dx2dy2 + 5*dx2dz2 + 16*dy4 - 10*dy2dz2 + 2*dz4; - double b5 = 20*dx4 - 15*dx2dy2 - 35*dx2dz2 - 40*dy4 + 115*dy2dz2 - 32*dz4; - double b6 = 10*dx4 + 15*dx2dy2 - 40*dx2dz2 + 10*dy4 - 40*dy2dz2 + 32*dz4; - term7 = b1*sx4 + b2*sx2*sy2 + b3*sx2*sz2 - + b4*sy4 + b5*sy2*sz2 + b6*sz4; - term7 *= (7.0/16.0); + double b1 = 16 * dx4 - 30 * dx2dy2 - 10 * dx2dz2 + 10 * dy4 + 5 * dy2dz2 + 2 * dz4; + double b2 = -40 * dx4 + 105 * dx2dy2 - 5 * dx2dz2 - 40 * dy4 - 5 * dy2dz2 + 4 * dz4; + double b3 = -40 * dx4 - 15 * dx2dy2 + 115 * dx2dz2 + 20 * dy4 - 35 * dy2dz2 - 32 * dz4; + double b4 = 10 * dx4 - 30 * dx2dy2 + 5 * dx2dz2 + 16 * dy4 - 10 * dy2dz2 + 2 * dz4; + double b5 = 20 * dx4 - 15 * dx2dy2 - 35 * dx2dz2 - 40 * dy4 + 115 * dy2dz2 - 32 * dz4; + double b6 = 10 * dx4 + 15 * dx2dy2 - 40 * dx2dz2 + 10 * dy4 - 40 * dy2dz2 + 32 * dz4; + term7 = b1 * sx4 + b2 * sx2 * sy2 + b3 * sx2 * sz2 + b4 * sy4 + b5 * sy2 * sz2 + b6 * sz4; + 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 * WIDE_PI * R2)) * (((term7 / R2 + term5) / R2 + term3) / (R2 * R)); // Error should be of order 1/R^9 return Nxy; } - double DemagNyyAsymptotic(double x, double y, double z, - double dx,double dy,double dz) -{ - return DemagNxxAsymptotic(y,x,z,dy,dx,dz); + double dx, double dy, double dz) { + return DemagNxxAsymptotic(y, x, z, dy, dx, dz); } double DemagNzzAsymptotic(double x, double y, double z, - double dx,double dy,double dz) -{ - return DemagNxxAsymptotic(z,y,x,dz,dy,dx); + double dx, double dy, double dz) { + return DemagNxxAsymptotic(z, y, x, dz, dy, dx); } double DemagNxzAsymptotic(double x, double y, double z, - double dx,double dy,double dz) -{ - return DemagNxyAsymptotic(x,z,y,dx,dz,dy); + double dx, double dy, double dz) { + return DemagNxyAsymptotic(x, z, y, dx, dz, dy); } double DemagNyzAsymptotic(double x, double y, double z, - double dx,double dy,double dz) -{ - return DemagNxyAsymptotic(y,z,x,dy,dz,dx); + double dx, double dy, double dz) { + return DemagNxyAsymptotic(y, z, x, dy, dz, dx); } diff --git a/fidimag/common/dipolar/demagcoef.h b/fidimag/common/dipolar/demagcoef.h index 94e34483..8729cce5 100644 --- a/fidimag/common/dipolar/demagcoef.h +++ b/fidimag/common/dipolar/demagcoef.h @@ -1,9 +1,9 @@ /* 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 @@ -39,10 +39,6 @@ to copyright protection within the United States. */ - - - - /* FILE: demagcoef.h -*-Mode: c++-*- * * Demag coefficients. @@ -73,83 +69,75 @@ to copyright protection within the United States. * */ +#include -double Newell_f(double x,double y,double z); -double Newell_g(double x,double y,double z); +double Newell_f(double x, double y, double z); +double Newell_g(double x, double y, double z); double CalculateSDA00(double x, double y, double z, - double dx,double dy,double dz); + double dx, double dy, double dz); inline double CalculateSDA11(double x, double y, double z, - double dx,double dy,double dz) -{ return CalculateSDA00(y,x,z,dy,dx,dz); } + double dx, double dy, double dz) { return CalculateSDA00(y, x, z, dy, dx, dz); } inline double CalculateSDA22(double x, double y, double z, - double dx,double dy,double dz) -{ return CalculateSDA00(z,y,x,dz,dy,dx); } - + double dx, double dy, double dz) { return CalculateSDA00(z, y, x, dz, dy, dx); } double CalculateSDA01(double x, double y, double z, - double dx,double dy,double dz); + double dx, double dy, double dz); inline double CalculateSDA02(double x, double y, double z, - double dx,double dy,double dz) -{ return CalculateSDA01(x,z,y,dx,dz,dy); } + double dx, double dy, double dz) { return CalculateSDA01(x, z, y, dx, dz, dy); } inline double CalculateSDA12(double x, double y, double z, - double dx,double dy,double dz) -{ return CalculateSDA01(y,z,x,dy,dz,dx); } + double dx, double dy, double dz) { return CalculateSDA01(y, z, x, dy, dz, dx); } - -double SelfDemagNx(double xsize,double ysize,double zsize); -double SelfDemagNy(double xsize,double ysize,double zsize); -double SelfDemagNz(double xsize,double ysize,double zsize); +double SelfDemagNx(double xsize, double ysize, double zsize); +double SelfDemagNy(double xsize, double ysize, double zsize); +double SelfDemagNz(double xsize, double ysize, double zsize); double DemagNxxAsymptotic(double x, double y, double z, - double dx,double dy,double dz); + double dx, double dy, double dz); double DemagNxyAsymptotic(double x, double y, double z, - double dx,double dy,double dz); + double dx, double dy, double dz); double DemagNyyAsymptotic(double x, double y, double z, - double dx,double dy,double dz); + double dx, double dy, double dz); double DemagNzzAsymptotic(double x, double y, double z, - double dx,double dy,double dz); + double dx, double dy, double dz); double DemagNxzAsymptotic(double x, double y, double z, - double dx,double dy,double dz); + double dx, double dy, double dz); double DemagNyzAsymptotic(double x, double y, double z, - double dx,double dy,double dz); - - + double dx, double dy, double dz); //////////////////////////////////////////////////////////////////////////// // Routines to do accurate summation -static int AS_Compare(const void* px,const void* py) -{ +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; + double x = fabs(*((const double *)px)); + double y = fabs(*((const double *)py)); + if (x < y) + return 1; + if (x > y) + return -1; return 0; } - -static double -AccurateSum(int n,double *arr) -{ +static double +AccurateSum(int n, double *arr) { // Order by decreasing magnitude - qsort(arr,n,sizeof(double),AS_Compare); + qsort(arr, n, sizeof(double), AS_Compare); // Add up using doubly compensated summation. If necessary, mark // variables these "volatile" to protect against problems arising @@ -157,22 +145,22 @@ AccurateSum(int n,double *arr) // 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 +#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; } - +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 + Tensor_xx, + Tensor_yy, + Tensor_zz, + Tensor_xy, + Tensor_xz, + Tensor_yz }; //========================================== -//used for demag +// 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; + 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); + 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); @@ -82,28 +86,26 @@ void compute_fields(fft_demag_plan *restrict plan, double *restrict spin, double 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 +// used for sode typedef struct { - int nxyz; + int nxyz; - double dt; - double T; - double gamma; - double *mu_s; - double coeff; - double Q; + double dt; + double T; + double gamma; + double *mu_s; + double coeff; + double Q; - double theta; - double theta1; - double theta2; + double theta; + double theta1; + double theta2; - double *dm1; - double *dm2; - double *eta; + double *dm1; + double *dm2; + double *eta; } ode_solver; @@ -111,7 +113,6 @@ void init_solver(ode_solver *s, double k_B, double theta, int nxyz, double dt, d 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); + 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); - + double *alpha, double *mu_s_inv, int *pins); diff --git a/fidimag/common/dipolar/tensor_2dpbc.h b/fidimag/common/dipolar/tensor_2dpbc.h index e2d41ae5..bbb64f4d 100644 --- a/fidimag/common/dipolar/tensor_2dpbc.h +++ b/fidimag/common/dipolar/tensor_2dpbc.h @@ -1,128 +1,113 @@ //------------------------------------------------------------------------------ -inline double RR(double x,double y,double z) -{return x*x+y*y+z*z;} +inline double RR(double x, double y, double z) { return x * x + y * y + z * z; } -inline double R(double x,double y,double z) -{return sqrt(x*x+y*y+z*z);} +inline double R(double x, double y, double z) { return sqrt(x * x + y * y + z * z); } -//------------------------------------------------------------------------------ +//------------------------------------------------------------------------------ //------------------------------------------------------------------------------ -//about Nxxinf - -inline double -fxx(double x,double y,double z) -{ - double t1,t2,t3,t5; - t1 = x*x; - t2 = y*y; - t3 = z*z; - t5 = sqrt(t1 + t2 + t3); - return y * x / (t1 + t3) / t5; +// about Nxxinf + +inline double +fxx(double x, double y, double z) { + double t1, t2, t3, t5; + t1 = x * x; + t2 = y * y; + t3 = z * z; + t5 = sqrt(t1 + t2 + t3); + return y * x / (t1 + t3) / t5; } - -inline double Nxxinf(double x,double y,double z,double X0,double Y0) { - return fxx(x+X0,y-Y0,z)+fxx(x-X0,y+Y0,z)-fxx(x+X0,y+Y0,z)-fxx(x-X0,y-Y0,z); +inline double Nxxinf(double x, double y, double z, double X0, double Y0) { + return fxx(x + X0, y - Y0, z) + fxx(x - X0, y + Y0, z) - fxx(x + X0, y + Y0, z) - fxx(x - X0, y - Y0, z); } //------------------------------------------------------------------------------ -//about Nyyinf +// about Nyyinf -inline double Nyyinf(double x,double y,double z,double X0,double Y0){ - return Nxxinf(y,x,z,Y0,X0); +inline double Nyyinf(double x, double y, double z, double X0, double Y0) { + return Nxxinf(y, x, z, Y0, X0); } - //------------------------------------------------------------------------------ -//about Nzzinf - -inline double fzz(double x,double y,double z){ - double t1,t2,t3,t5; - t1 = x*x; - t2 = y*y; - t3 = z*z; - t5 = sqrt(t1 + t2 + t3); - return y*x*(t1+t2+2*t3)/((t1 + t3)*(t2+t3)*t5); +// about Nzzinf + +inline double fzz(double x, double y, double z) { + double t1, t2, t3, t5; + t1 = x * x; + t2 = y * y; + t3 = z * z; + t5 = sqrt(t1 + t2 + t3); + return y * x * (t1 + t2 + 2 * t3) / ((t1 + t3) * (t2 + t3) * t5); } - -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); +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); } - - //------------------------------------------------------------------------------ -//about Nxyinf -inline double fxy(double x,double y,double z){ - return 1.0/R(x,y,z); +// about Nxyinf +inline double fxy(double x, double y, double z) { + return 1.0 / R(x, y, z); } - -inline double Nxyinf(double x,double y,double z,double X0,double Y0){ - return (fxy(x-X0,y-Y0,z)+fxy(x+X0,y+Y0,z)-fxy(x-X0,y+Y0,z)-fxy(x+X0,y-Y0,z)); +inline double Nxyinf(double x, double y, double z, double X0, double Y0) { + return (fxy(x - X0, y - Y0, z) + fxy(x + X0, y + Y0, z) - fxy(x - X0, y + Y0, z) - fxy(x + X0, y - Y0, z)); } - //------------------------------------------------------------------------------ -//about Nxzinf -inline double fxz(double x,double y,double z){ - double t2,t3,t4,t6; - t2=x*x; - t3=y*y; - t4=z*z; - t6 = sqrt(t2 + t3 + t4); - return y * z / (t6 * (t2 + t4)); +// about Nxzinf +inline double fxz(double x, double y, double z) { + double t2, t3, t4, t6; + t2 = x * x; + t3 = y * y; + t4 = z * z; + t6 = sqrt(t2 + t3 + t4); + return y * z / (t6 * (t2 + t4)); } - -inline double Nxzinf(double x,double y,double z,double X0,double Y0){ - return fxz(x+X0,y-Y0,z)+fxz(x-X0,y+Y0,z)-fxz(x+X0,y+Y0,z)-fxz(x-X0,y-Y0,z); +inline double Nxzinf(double x, double y, double z, double X0, double Y0) { + return fxz(x + X0, y - Y0, z) + fxz(x - X0, y + Y0, z) - fxz(x + X0, y + Y0, z) - fxz(x - X0, y - Y0, z); } //------------------------------------------------------------------------------ -//about Nyzinf +// about Nyzinf -inline double Nyzinf(double x,double y,double z,double X0,double Y0){ - return Nxzinf(y,x,z,Y0,X0); +inline double Nyzinf(double x, double y, double z, double X0, double Y0) { + return Nxzinf(y, x, z, Y0, X0); } - //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ -//about Nxxdipole -inline double -DemagNxxDipolar(double x,double y,double z) -{ -double t1 = x*x; -double t3 = y*y; -double t4 = z*z; -double t6 = t1 + t3 + t4;; -double t7 = t6*t6; -double t8 = sqrt(t6); -return -(2.0 * t1 - t3 - t4) / (t8 * t7); +// about Nxxdipole +inline double +DemagNxxDipolar(double x, double y, double z) { + double t1 = x * x; + double t3 = y * y; + double t4 = z * z; + double t6 = t1 + t3 + t4; + ; + double t7 = t6 * t6; + double t8 = sqrt(t6); + return -(2.0 * t1 - t3 - t4) / (t8 * t7); } -inline double -DemagNxyDipolar(double x,double y,double z) -{ -double t6 = RR(x,y,z); -double t7 = t6*t6; -double t8 = sqrt(t6); -return -3.0*x*y/(t8*t7); +inline double +DemagNxyDipolar(double x, double y, double z) { + double t6 = RR(x, y, z); + double t7 = t6 * t6; + double t8 = sqrt(t6); + return -3.0 * x * y / (t8 * t7); } //----------------------------------------------------------------------------- -double DemagTensorNormal(enum Type_Nij comp,double x,double y,double z,double a,double b,double c); - -double DemagTensorAsymptotic(enum Type_Nij comp,double x,double y,double z,double a,double b,double c); +double DemagTensorNormal(enum Type_Nij comp, double x, double y, double z, double a, double b, double c); -double DemagTensorDipolar(enum Type_Nij comp,double x,double y,double z); +double DemagTensorAsymptotic(enum Type_Nij comp, double x, double y, double z, double a, double b, double c); -double DemagTensorInfinite(enum Type_Nij comp,double x,double y,double z,double X0,double Y0); +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/fidimag/common/lib/common_clib.h b/fidimag/common/lib/common_clib.h index d7a024ab..5d0cb9c7 100644 --- a/fidimag/common/lib/common_clib.h +++ b/fidimag/common/lib/common_clib.h @@ -2,7 +2,7 @@ #define __CLIB__ #include -//#include +// #include #define WIDE_PI 3.1415926535897932384626433832795L // ---------------------------------------------------------------------------- @@ -10,35 +10,36 @@ /* 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; + 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; + 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; + 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; +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; + if (pins[id] > 0) + continue; mm = sqrt(m[i] * m[i] + m[j] * m[j] + m[k] * m[k]); - if(mm > 0) { + if (mm > 0) { mm = 1 / mm; m[i] *= mm; m[j] *= mm; m[k] *= mm; } - } + } } // ---------------------------------------------------------------------------- @@ -67,13 +68,12 @@ void llg_stt_cpp(double *restrict dm_dt, double *restrict m, double *restrict h, // ---------------------------------------------------------------------------- // 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); +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/llg.c b/fidimag/common/lib/llg.c index ce438fdf..357b04f6 100644 --- a/fidimag/common/lib/llg.c +++ b/fidimag/common/lib/llg.c @@ -18,46 +18,46 @@ * m is the magnetisation vector and H_eff is the effective field. * * In our calculation, we usually compute: - + * m x (m x H_eff) = ( m * H_eff) m - (m * m) H_eff - + * Then, if we define the perpendicular component of the effective * field as: - + * H_perp = H_eff - ( m * H_eff) m - + * we can write - + * dm -gamma * ---- = -------- ( m X H_perp - a * H_perp ) * dt 2 * ( 1 + a ) - + * since m X m = 0 and, theoretically, m * m = 1 (at zero temperature m has * fixed length). However, for the second term to the right hand side, * proportional to a, it is better to keep the m X ( m x H_eff ) term with (m * * m) for stability purposes, thus we use * - + * H_perp = (m * m) H_eff - ( m * H_eff) m - - * for both terms. + + * for both terms. * * Additionally, to preserve the magnetisation length, we need to correct the * dm / dt term for every time step, adding the following term to the right * hand size expression of the LLG: - + * dm dm 2 * ---- ---> ---- + c * ( 1 - m ) m - * dt dt - + * dt dt + * with _________ * / 2 - * c = 6 * / ( dm ) + * c = 6 * / ( dm ) * / ( ---- ) * \/ ( dt ) - + * The correction must be introduced since numerically, m can change length * during a step of the integration (which would not occur if the integration * step is infinitely small), deviating from the real answer. If we just @@ -76,18 +76,18 @@ * dm/dt changes the time scaling by a factor proportional to 1/t, therefore in * the future we could try to estimate its * influence with more mathematical/numerical rigour and analyse an optimal - * value for the prefactor (6). + * value for the prefactor (6). */ void llg_rhs(double *restrict dm_dt, double *restrict m, double *restrict h, double *restrict alpha, int *restrict pins, - double gamma, int n, int do_precession, double default_c) { + double gamma, int n, int do_precession, double default_c) { 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) +#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 @@ -109,7 +109,7 @@ void llg_rhs(double *restrict dm_dt, double *restrict m, double *restrict h, dou mm = m[i] * m[i] + m[j] * m[j] + m[k] * m[k]; mh = m[i] * h[i] + m[j] * h[j] + m[k] * h[k]; - // Usually, m is normalised, i.e., mm=1; + // 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 @@ -129,7 +129,7 @@ void llg_rhs(double *restrict dm_dt, double *restrict m, double *restrict h, dou double mth0 = 0, mth1 = 0, mth2 = 0; // The first term: m x H_eff = m x H_perp - if (do_precession){ + if (do_precession) { 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); @@ -146,70 +146,61 @@ void llg_rhs(double *restrict dm_dt, double *restrict m, double *restrict h, dou // (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){ + 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] - ); + 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); + // printf("%0.15g %0.15g\n", c, default_c); // Correct the RHS term to keep m normalised dm_dt[i] += c * (1 - mm) * m[i]; dm_dt[j] += c * (1 - mm) * m[j]; dm_dt[k] += c * (1 - mm) * m[k]; - } - } - 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) { + 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) + // #pragma omp parallel for private(i,j,k,coeff,mm, mh, c, hpi,hpj,hpk) for (int id = 0; id < n; id++) { int i = 3 * id; int j = i + 1; int k = i + 2; - - if (pins[i]>0){ - continue; + if (pins[i] > 0) { + continue; } - double coeff = -gamma/(1.0+alpha[i]*alpha[i]); + double coeff = -gamma / (1.0 + alpha[i] * alpha[i]); - if (do_precession){ - jtn[i] = coeff*(cross_x(mp[i],mp[j],mp[k],h[i],h[j],h[k])+cross_x(m[i],m[j],m[k],hp[i],hp[j],hp[k])); - jtn[j] = coeff*(cross_y(mp[i],mp[j],mp[k],h[i],h[j],h[k])+cross_y(m[i],m[j],m[k],hp[i],hp[j],hp[k])); - jtn[k] = coeff*(cross_z(mp[i],mp[j],mp[k],h[i],h[j],h[k])+cross_z(m[i],m[j],m[k],hp[i],hp[j],hp[k])); - }else{ + if (do_precession) { + jtn[i] = coeff * (cross_x(mp[i], mp[j], mp[k], h[i], h[j], h[k]) + cross_x(m[i], m[j], m[k], hp[i], hp[j], hp[k])); + jtn[j] = coeff * (cross_y(mp[i], mp[j], mp[k], h[i], h[j], h[k]) + cross_y(m[i], m[j], m[k], hp[i], hp[j], hp[k])); + jtn[k] = coeff * (cross_z(mp[i], mp[j], mp[k], h[i], h[j], h[k]) + cross_z(m[i], m[j], m[k], hp[i], hp[j], hp[k])); + } else { jtn[i] = 0; jtn[j] = 0; jtn[k] = 0; } - double mm = m[i]*m[i] + m[j]*m[j] + m[k]*m[k]; - double mh = m[i]*h[i] + m[j]*h[j] + m[k]*h[k]; - double mhp = m[i]*hp[i] + m[j]*hp[j] + m[k]*hp[k]; - double mph = mp[i]*h[i] + mp[j]*h[j] + mp[k]*h[k]; - double mmp = m[i]*mp[i] + m[j]*mp[j] + m[k]*mp[k]; - - jtn[i] += alpha[i]*coeff*((mph+mhp)*m[i]+mh*mp[i]-2*mmp*h[i]-mm*hp[i]); - jtn[j] += alpha[i]*coeff*((mph+mhp)*m[j]+mh*mp[j]-2*mmp*h[j]-mm*hp[j]); - jtn[k] += alpha[i]*coeff*((mph+mhp)*m[k]+mh*mp[k]-2*mmp*h[k]-mm*hp[k]); - - - if (default_c>0){ - jtn[i] += default_c *((1-mm)*mp[i]-2*mmp*m[i]); - jtn[j] += default_c *((1-mm)*mp[j]-2*mmp*m[j]); - jtn[k] += default_c *((1-mm)*mp[k]-2*mmp*m[k]); - } + double mm = m[i] * m[i] + m[j] * m[j] + m[k] * m[k]; + double mh = m[i] * h[i] + m[j] * h[j] + m[k] * h[k]; + double mhp = m[i] * hp[i] + m[j] * hp[j] + m[k] * hp[k]; + double mph = mp[i] * h[i] + mp[j] * h[j] + mp[k] * h[k]; + double mmp = m[i] * mp[i] + m[j] * mp[j] + m[k] * mp[k]; - } + jtn[i] += alpha[i] * coeff * ((mph + mhp) * m[i] + mh * mp[i] - 2 * mmp * h[i] - mm * hp[i]); + jtn[j] += alpha[i] * coeff * ((mph + mhp) * m[j] + mh * mp[j] - 2 * mmp * h[j] - mm * hp[j]); + jtn[k] += alpha[i] * coeff * ((mph + mhp) * m[k] + mh * mp[k] - 2 * mmp * h[k] - mm * hp[k]); + if (default_c > 0) { + jtn[i] += default_c * ((1 - mm) * mp[i] - 2 * mmp * m[i]); + jtn[j] += default_c * ((1 - mm) * mp[j] - 2 * mmp * m[j]); + jtn[k] += default_c * ((1 - mm) * mp[k] - 2 * mmp * m[k]); + } + } } - diff --git a/fidimag/common/lib/steepest_descent.c b/fidimag/common/lib/steepest_descent.c index f1d860b0..fb4ef075 100644 --- a/fidimag/common/lib/steepest_descent.c +++ b/fidimag/common/lib/steepest_descent.c @@ -1,12 +1,12 @@ #include "common_clib.h" #include "math.h" -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_update_spin(double *spin, double *spin_last, double *magnetisation, + double *mxH, double *mxmxH, double *mxmxH_last, double tau, + int *pins, int n) { - // Update the field ------------------------------------------------------- - #pragma omp parallel for +// Update the field ------------------------------------------------------- +#pragma omp parallel for for (int i = 0; i < n; i++) { if (magnetisation[i] > 0.0) { @@ -33,26 +33,25 @@ void sd_update_spin (double *spin, double *spin_last, double *magnetisation, factor_minus = 4 - tau * tau * mxH_sq; for (int j = 0; j < 3; j++) { - new_spin[j] = factor_minus * spin[spin_idx + j] - - 4 * tau * mxmxH[spin_idx + j]; + new_spin[j] = factor_minus * spin[spin_idx + j] - 4 * tau * mxmxH[spin_idx + j]; new_spin[j] = new_spin[j] / factor_plus; spin[spin_idx + j] = new_spin[j]; } - } // close if Ms or mu_s > 0 - } // close for + } // close if Ms or mu_s > 0 + } // close for normalise(spin, pins, 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) { +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) { double res; double num = 0; double den = 0; double sign; - #pragma omp parallel for reduction(+: num,den) +#pragma omp parallel for reduction(+ : num, den) for (int i = 0; i < n; i++) { if (magnetisation[i] > 0.0) { int spin_idx; @@ -67,10 +66,10 @@ void sd_compute_step (double *spin, double *spin_last, double *magnetisation, do // Compute the torques - mxH[spin_idx] = cross_x(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], - field[spin_idx], - field[spin_idx + 1], - field[spin_idx + 2]); + mxH[spin_idx] = cross_x(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], + field[spin_idx], + field[spin_idx + 1], + field[spin_idx + 2]); mxH[spin_idx + 1] = cross_y(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], field[spin_idx], @@ -82,8 +81,8 @@ void sd_compute_step (double *spin, double *spin_last, double *magnetisation, do field[spin_idx + 1], field[spin_idx + 2]); - mxmxH[spin_idx] = cross_x(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], - mxH[spin_idx], mxH[spin_idx + 1], mxH[spin_idx + 2]); + mxmxH[spin_idx] = cross_x(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], + mxH[spin_idx], mxH[spin_idx + 1], mxH[spin_idx + 2]); mxmxH[spin_idx + 1] = cross_y(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], mxH[spin_idx], mxH[spin_idx + 1], mxH[spin_idx + 2]); mxmxH[spin_idx + 2] = cross_z(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], @@ -100,23 +99,21 @@ void sd_compute_step (double *spin, double *spin_last, double *magnetisation, do num += ds[j] * ds[j]; den += ds[j] * dy[j]; } - } - else { + } else { for (int j = 0; j < 3; j++) { num += ds[j] * dy[j]; den += dy[j] * dy[j]; } } - } // Close if Ms or mu_s > 0 - } // close for + } // Close if Ms or mu_s > 0 + } // close for // Criteria for the evaluation of tau is in line 96 of: - //https://github.com/MicroMagnum/MicroMagnum/blob/minimizer/src/magnum/micromagnetics/micro_magnetics_solver.py + // https://github.com/MicroMagnum/MicroMagnum/blob/minimizer/src/magnum/micromagnetics/micro_magnetics_solver.py if (den == 0.0) { res = tmax; - } - else { + } else { res = num / den; } diff --git a/fidimag/common/lib/stt.c b/fidimag/common/lib/stt.c index b2faeb10..9b0c979a 100644 --- a/fidimag/common/lib/stt.c +++ b/fidimag/common/lib/stt.c @@ -52,16 +52,16 @@ * */ 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) { + double dx, double dy, double dz, int *restrict ngbs, int n) { - //#pragma omp parallel for - for (int i = 0; i < 3 * n; i++) { - field[i] = 0; - } + // #pragma omp parallel for + for (int i = 0; i < 3 * n; i++) { + field[i] = 0; + } - #pragma omp parallel for +#pragma omp parallel for /* Iterate through every lattice site */ - for (int i = 0; i < n; i++){ + for (int i = 0; i < n; i++) { /* Starting index for the NNs of the i-th site * i+0, i+1, i+2, i+3 ... --> -x, +x, -y, +y ... @@ -76,17 +76,17 @@ void compute_stt_field_c(double *restrict spin, double *restrict field, double * * In the latest case, make factor_x equal to zero to avoid * summing field to the for loop */ - if(ngbs[nn_i] != -1 && ngbs[nn_i + 1] != -1) { + if (ngbs[nn_i] != -1 && ngbs[nn_i + 1] != -1) { factor_x = 2; nn_x1 = ngbs[nn_i]; nn_x2 = ngbs[nn_i + 1]; // Here there is no NN to the right so we make f(x) - f(x-1) - } else if(ngbs[nn_i + 1] == -1 && ngbs[nn_i] != -1){ + } else if (ngbs[nn_i + 1] == -1 && ngbs[nn_i] != -1) { factor_x = 1; nn_x1 = ngbs[nn_i]; nn_x2 = i; // Here there is no NN to the left so we make f(x + 1) - f(x) - } else if(ngbs[nn_i] == -1 && ngbs[nn_i + 1] != -1){ + } else if (ngbs[nn_i] == -1 && ngbs[nn_i + 1] != -1) { factor_x = 1; nn_x1 = i; nn_x2 = ngbs[nn_i + 1]; @@ -99,23 +99,22 @@ void compute_stt_field_c(double *restrict spin, double *restrict field, double * * jx is a scalar field, so it only has n entries * This calculation is: jx[i] * d m[i] / dx * */ - if (factor_x){ - for(int j = 0; j < 3; j++){ - field[3 * i + j] += jx[i] * (spin[3 * nn_x2 + j] - - spin[3 * nn_x1 + j]) / (factor_x * dx); + if (factor_x) { + for (int j = 0; j < 3; j++) { + field[3 * i + j] += jx[i] * (spin[3 * nn_x2 + j] - spin[3 * nn_x1 + j]) / (factor_x * dx); } } // We do the same along the y direction - if(ngbs[nn_i + 2] != -1 && ngbs[nn_i + 3] != -1) { + if (ngbs[nn_i + 2] != -1 && ngbs[nn_i + 3] != -1) { factor_y = 2; nn_y1 = ngbs[nn_i + 2]; nn_y2 = ngbs[nn_i + 3]; - } else if(ngbs[nn_i + 3] == -1 && ngbs[nn_i + 2] != -1){ + } else if (ngbs[nn_i + 3] == -1 && ngbs[nn_i + 2] != -1) { factor_y = 1; nn_y1 = ngbs[nn_i + 2]; nn_y2 = i; - } else if(ngbs[nn_i + 2] == -1 && ngbs[nn_i + 3] != -1){ + } else if (ngbs[nn_i + 2] == -1 && ngbs[nn_i + 3] != -1) { factor_y = 1; nn_y1 = i; nn_y2 = ngbs[nn_i + 3]; @@ -123,24 +122,22 @@ void compute_stt_field_c(double *restrict spin, double *restrict field, double * factor_y = 0; } - if (factor_y){ - for(int j = 0; j < 3; j++){ - field[3 * i + j] += jy[i] * (spin[3 * nn_y2 + j] - - spin[3 * nn_y1 + j]) / (factor_y * dy); + if (factor_y) { + for (int j = 0; j < 3; j++) { + field[3 * i + j] += jy[i] * (spin[3 * nn_y2 + j] - spin[3 * nn_y1 + j]) / (factor_y * dy); } } - // We do the same along the z direction - if(ngbs[nn_i + 4] >= 0 && ngbs[nn_i + 5] >= 0) { + if (ngbs[nn_i + 4] >= 0 && ngbs[nn_i + 5] >= 0) { factor_z = 2; nn_z1 = ngbs[nn_i + 4]; nn_z2 = ngbs[nn_i + 5]; - } else if(ngbs[nn_i + 4] >= 0 && ngbs[nn_i + 5] < 0){ + } else if (ngbs[nn_i + 4] >= 0 && ngbs[nn_i + 5] < 0) { factor_z = 1; nn_z1 = ngbs[nn_i + 4]; nn_z2 = i; - } else if(ngbs[nn_i + 4] < 0 && ngbs[nn_i + 5] >= 0 ){ + } else if (ngbs[nn_i + 4] < 0 && ngbs[nn_i + 5] >= 0) { factor_z = 1; nn_z1 = i; nn_z2 = ngbs[nn_i + 5]; @@ -148,142 +145,127 @@ void compute_stt_field_c(double *restrict spin, double *restrict field, double * factor_z = 0; } - if (factor_z){ - for(int j = 0; j < 3; j++){ - field[3 * i + j] += jz[i] * (spin[3 * nn_z2 + j] - - spin[3 * nn_z1 + j]) / (factor_z * dz); + if (factor_z) { + for (int j = 0; j < 3; j++) { + field[3 * i + j] += jz[i] * (spin[3 * nn_z2 + j] - spin[3 * nn_z1 + j]) / (factor_z * dz); } } - } } - 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) { + double *restrict alpha, double beta, double u0, double gamma, int n) { - #pragma omp parallel for - for (int index = 0; index < n; index++) { - int i = 3 * index; - int j = 3 * index + 1; - int k = 3 * index + 2; +#pragma omp parallel for + for (int index = 0; index < n; index++) { + int i = 3 * index; + int j = 3 * index + 1; + int k = 3 * index + 2; - double coeff = -gamma / (1 + alpha[index] * alpha[index]); + double coeff = -gamma / (1 + alpha[index] * alpha[index]); - double mm = m[i] * m[i] + m[j] * m[j] + m[k] * m[k]; - double mh = m[i] * h[i] + m[j] * h[j] + m[k] * h[k]; + double mm = m[i] * m[i] + m[j] * m[j] + m[k] * m[k]; + double mh = m[i] * h[i] + m[j] * h[j] + m[k] * h[k]; - //hp=mm.h-mh.m=-mx(mxh) - double hpi = mm*h[i] - mh*m[i]; - double hpj = mm*h[j] - mh*m[j]; - double hpk = mm*h[k] - mh*m[k]; + // hp=mm.h-mh.m=-mx(mxh) + double hpi = mm * h[i] - mh * m[i]; + double hpj = mm * h[j] - mh * m[j]; + double hpk = mm * h[k] - mh * m[k]; - double mth0 = cross_x(m[i], m[j], m[k], hpi, hpj, hpk); - double mth1 = cross_y(m[i], m[j], m[k], hpi, hpj, hpk); - double mth2 = cross_z(m[i], m[j], m[k], hpi, hpj, hpk); + double mth0 = cross_x(m[i], m[j], m[k], hpi, hpj, hpk); + double mth1 = cross_y(m[i], m[j], m[k], hpi, hpj, hpk); + double mth2 = cross_z(m[i], m[j], m[k], hpi, hpj, hpk); - dm_dt[i] = coeff * (mth0 - hpi * alpha[index]); - dm_dt[j] = coeff * (mth1 - hpj * alpha[index]); - dm_dt[k] = coeff * (mth2 - hpk * alpha[index]); + dm_dt[i] = coeff * (mth0 - hpi * alpha[index]); + dm_dt[j] = coeff * (mth1 - hpj * alpha[index]); + dm_dt[k] = coeff * (mth2 - hpk * alpha[index]); - //the above part is standard LLG equation. + // the above part is standard LLG equation. - double coeff_stt = u0 / (1 + alpha[index] * alpha[index]); + double coeff_stt = u0 / (1 + alpha[index] * alpha[index]); - double mht = m[i] * h_stt[i] + m[j] * h_stt[j] + m[k] * h_stt[k]; + double mht = m[i] * h_stt[i] + m[j] * h_stt[j] + m[k] * h_stt[k]; - hpi = mm*h_stt[i] - mht * m[i]; - hpj = mm*h_stt[j] - mht * m[j]; - hpk = mm*h_stt[k] - mht * m[k]; + hpi = mm * h_stt[i] - mht * m[i]; + hpj = mm * h_stt[j] - mht * m[j]; + hpk = mm * h_stt[k] - mht * m[k]; - 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); + 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); - dm_dt[i] += coeff_stt * ((1 + alpha[index] * beta) * hpi - - (beta - alpha[index]) * mth0); - dm_dt[j] += coeff_stt * ((1 + alpha[index] * beta) * hpj - - (beta - alpha[index]) * mth1); - dm_dt[k] += coeff_stt * ((1 + alpha[index] * beta) * hpk - - (beta - alpha[index]) * mth2); + dm_dt[i] += coeff_stt * ((1 + alpha[index] * beta) * hpi - (beta - alpha[index]) * mth0); + dm_dt[j] += coeff_stt * ((1 + alpha[index] * beta) * hpj - (beta - alpha[index]) * mth1); + dm_dt[k] += coeff_stt * ((1 + alpha[index] * beta) * hpk - (beta - alpha[index]) * mth2); - double c = 6 * sqrt(dm_dt[i] * dm_dt[i] + - dm_dt[j] * dm_dt[j] + - dm_dt[k]* dm_dt[k]); - - dm_dt[i] += c * (1 - mm) * m[i]; - dm_dt[j] += c * (1 - mm) * m[j]; - dm_dt[k] += c * (1 - mm) * m[k]; - - } + double c = 6 * sqrt(dm_dt[i] * dm_dt[i] + + dm_dt[j] * dm_dt[j] + + dm_dt[k] * dm_dt[k]); + dm_dt[i] += c * (1 - mm) * m[i]; + dm_dt[j] += c * (1 - mm) * m[j]; + dm_dt[k] += c * (1 - mm) * m[k]; + } } - 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) { - - #pragma omp parallel for - for (int index = 0; index < n; index++) { - int i = 3 * index; - int j = 3 * index + 1; - int k = 3 * index + 2; - - if (pins[index]>0){ - dm_dt[i] = 0; - dm_dt[j] = 0; - dm_dt[k] = 0; - continue; - } - - double coeff = -gamma / (1 + alpha[index] * alpha[index]); - - double mm = m[i] * m[i] + m[j] * m[j] + m[k] * m[k]; - double mh = m[i] * h[i] + m[j] * h[j] + m[k] * h[k]; - - //hp=mm.h-mh.m=-mx(mxh) - double hpi = mm*h[i] - mh*m[i]; - double hpj = mm*h[j] - mh*m[j]; - double hpk = mm*h[k] - mh*m[k]; + double *restrict alpha, int *restrict pins, double *restrict a_J, double beta, double gamma, int n) { + +#pragma omp parallel for + for (int index = 0; index < n; index++) { + int i = 3 * index; + int j = 3 * index + 1; + int k = 3 * index + 2; + + if (pins[index] > 0) { + dm_dt[i] = 0; + dm_dt[j] = 0; + dm_dt[k] = 0; + continue; + } - double mth0 = cross_x(m[i], m[j], m[k], hpi, hpj, hpk); - double mth1 = cross_y(m[i], m[j], m[k], hpi, hpj, hpk); - double mth2 = cross_z(m[i], m[j], m[k], hpi, hpj, hpk); + double coeff = -gamma / (1 + alpha[index] * alpha[index]); - dm_dt[i] = coeff * (mth0 - hpi * alpha[index]); - dm_dt[j] = coeff * (mth1 - hpj * alpha[index]); - dm_dt[k] = coeff * (mth2 - hpk * alpha[index]); + double mm = m[i] * m[i] + m[j] * m[j] + m[k] * m[k]; + double mh = m[i] * h[i] + m[j] * h[j] + m[k] * h[k]; - //the above part is standard LLG equation. + // hp=mm.h-mh.m=-mx(mxh) + double hpi = mm * h[i] - mh * m[i]; + double hpj = mm * h[j] - mh * m[j]; + double hpk = mm * h[k] - mh * m[k]; - double coeff_stt = a_J[index] / (1 + alpha[index] * alpha[index]); + double mth0 = cross_x(m[i], m[j], m[k], hpi, hpj, hpk); + double mth1 = cross_y(m[i], m[j], m[k], hpi, hpj, hpk); + double mth2 = cross_z(m[i], m[j], m[k], hpi, hpj, hpk); - double mp = m[i] * p[i] + m[j] * p[j] + m[k] * p[k]; + dm_dt[i] = coeff * (mth0 - hpi * alpha[index]); + dm_dt[j] = coeff * (mth1 - hpj * alpha[index]); + dm_dt[k] = coeff * (mth2 - hpk * alpha[index]); - hpi = mm*p[i] - mp * m[i]; - hpj = mm*p[j] - mp * m[j]; - hpk = mm*p[k] - mp * m[k]; + // the above part is standard LLG equation. - 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); + double coeff_stt = a_J[index] / (1 + alpha[index] * alpha[index]); + double mp = m[i] * p[i] + m[j] * p[j] + m[k] * p[k]; - dm_dt[i] += coeff_stt * ((1 + alpha[index] * beta) * hpi - - (beta - alpha[index]) * mth0); - dm_dt[j] += coeff_stt * ((1 + alpha[index] * beta) * hpj - - (beta - alpha[index]) * mth1); - dm_dt[k] += coeff_stt * ((1 + alpha[index] * beta) * hpk - - (beta - alpha[index]) * mth2); + hpi = mm * p[i] - mp * m[i]; + hpj = mm * p[j] - mp * m[j]; + hpk = mm * p[k] - mp * m[k]; - double c = 6 * sqrt(dm_dt[i] * dm_dt[i] + - dm_dt[j] * dm_dt[j] + - dm_dt[k]* dm_dt[k]); + 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); - dm_dt[i] += c * (1 - mm) * m[i]; - dm_dt[j] += c * (1 - mm) * m[j]; - dm_dt[k] += c * (1 - mm) * m[k]; + dm_dt[i] += coeff_stt * ((1 + alpha[index] * beta) * hpi - (beta - alpha[index]) * mth0); + dm_dt[j] += coeff_stt * ((1 + alpha[index] * beta) * hpj - (beta - alpha[index]) * mth1); + dm_dt[k] += coeff_stt * ((1 + alpha[index] * beta) * hpk - (beta - alpha[index]) * mth2); - } + double c = 6 * sqrt(dm_dt[i] * dm_dt[i] + + dm_dt[j] * dm_dt[j] + + dm_dt[k] * dm_dt[k]); + dm_dt[i] += c * (1 - mm) * m[i]; + dm_dt[j] += c * (1 - mm) * m[j]; + dm_dt[k] += c * (1 - mm) * m[k]; + } } diff --git a/fidimag/common/neb_method/nebm_geodesic_lib.c b/fidimag/common/neb_method/nebm_geodesic_lib.c index 7bc28736..bce5c6b5 100644 --- a/fidimag/common/neb_method/nebm_geodesic_lib.c +++ b/fidimag/common/neb_method/nebm_geodesic_lib.c @@ -1,12 +1,11 @@ #include "nebm_geodesic_lib.h" -#include "nebm_lib.h" #include "math.h" +#include "nebm_lib.h" double compute_geodesic_Vincenty(double *restrict A, double *restrict B, int n_dofs_image, int *restrict material, - int n_dofs_image_material - ) { + int n_dofs_image_material) { /* Compute the Geodesic distance between two images: A and B, of an energy * band. For this, we use Vincenty's formula, which is defined in @@ -24,14 +23,14 @@ double compute_geodesic_Vincenty(double *restrict A, double *restrict B, * the A or B arrays. Since we assume Cartesian * coordinates, the number of spins is just: * n_dofs_image / 3 - * + * * material :: An array of size (3 * n_spins), containing only 1 and 0s * For every spin, its value is repeated 3 times (the number * of dofs). The values of material is 1 where mu_s or M_s * are larger than zero. * This array is not necessary here since sites without * material do not change their magnetisation and thus they - * do not contribute to the Geodesic distance magnitude + * do not contribute to the Geodesic distance magnitude * * Vicenty's formula is defined by computing the distance between * corresponding spins of the images A and B. So, if we have m(A)_i as @@ -61,7 +60,7 @@ double compute_geodesic_Vincenty(double *restrict A, double *restrict B, // of A and B. The i-th spin components start at the 3 * i position // in the arrays // We do not sum sites without - for(int i = 0; i < n_spins; i++){ + for (int i = 0; i < n_spins; i++) { // We only need to know if the spin is in a site with material (Ms, mu_s > 0) // so we skip the material for m_y and m_z which are the same as // the material for m_x @@ -83,13 +82,12 @@ double compute_geodesic_Vincenty(double *restrict A, double *restrict B, } double compute_geodesic_GreatCircle(double *restrict A, double *restrict B, - int n_dofs_image, - int *restrict material, - int n_dofs_image_material - ) { + int n_dofs_image, + int *restrict material, + int n_dofs_image_material) { /* Compute the Geodesic distance between two images: A and B, of an energy - * band. + * band. * * INPUTS: * @@ -103,14 +101,14 @@ double compute_geodesic_GreatCircle(double *restrict A, double *restrict B, * the A or B arrays. Since we assume Cartesian * coordinates, the number of spins is just: * n_dofs_image / 3 - * + * * material :: An array of size (3 * n_spins), containing only 1 and 0s * For every spin, its value is repeated 3 times (the number * of dofs). The values of material is 1 where mu_s or M_s * are larger than zero. * This array is not necessary here since sites without * material do not change their magnetisation and thus they - * do not contribute to the Geodesic distance magnitude + * do not contribute to the Geodesic distance magnitude * */ @@ -124,7 +122,7 @@ double compute_geodesic_GreatCircle(double *restrict A, double *restrict B, // of A and B. The i-th spin components start at the 3 * i position // in the arrays // We do not sum sites without - for(int i = 0; i < n_spins; i++){ + for (int i = 0; i < n_spins; i++) { // We only need to know if the spin is in a site with material (Ms, mu_s > 0) // so we skip the material for m_y and m_z which are the same as // the material for m_x diff --git a/fidimag/common/neb_method/nebm_geodesic_lib.h b/fidimag/common/neb_method/nebm_geodesic_lib.h index 0de9c151..f0876e4b 100644 --- a/fidimag/common/neb_method/nebm_geodesic_lib.h +++ b/fidimag/common/neb_method/nebm_geodesic_lib.h @@ -4,11 +4,9 @@ double compute_geodesic_GreatCircle(double *restrict A, double *restrict B, int n_dofs_image, int *restrict material, - int n_dofs_image_material - ); + int n_dofs_image_material); double compute_geodesic_Vincenty(double *restrict A, double *restrict B, int n_dofs_image, int *restrict material, - int n_dofs_image_material - ); + int n_dofs_image_material); diff --git a/fidimag/common/neb_method/nebm_integrators.c b/fidimag/common/neb_method/nebm_integrators.c index efaa39d3..adc5fcce 100644 --- a/fidimag/common/neb_method/nebm_integrators.c +++ b/fidimag/common/neb_method/nebm_integrators.c @@ -1,7 +1,7 @@ #include "math.h" #include "nebm_lib.h" -#include #include +#include double step_Verlet_C(double *restrict forces, double *restrict forces_prev, @@ -13,51 +13,46 @@ double step_Verlet_C(double *restrict forces, double mass, int n_images, int n_dofs_image, - double (* update_field) (double, double *) - - ) { + double (*update_field)(double, double *) + +) { printf("Hellooooooo\n"); update_field(t, y); int im_idx; - for(int i = 1; i < n_images - 1; i++){ - + for (int i = 1; i < n_images - 1; i++) { + im_idx = n_dofs_image * i; - double * force = &forces[im_idx]; - double * force_prev = &forces_prev[im_idx]; - double * velocity = &velocities[im_idx]; - double * velocity_new = &velocities_new[im_idx]; + double *force = &forces[im_idx]; + double *force_prev = &forces_prev[im_idx]; + double *velocity = &velocities[im_idx]; + double *velocity_new = &velocities_new[im_idx]; double v_dot_f = 0; double f_dot_f = 0; - - for(int j = 0; i < n_dofs_image; j++){ - y[im_idx + j] += h * (velocity[j] - + (h / (2 * mass)) * force[j]); - velocity[j] = velocity_new[j] + - (h / (2 * mass)) * (force_prev[j] + force[j]); + for (int j = 0; i < n_dofs_image; j++) { + y[im_idx + j] += h * (velocity[j] + (h / (2 * mass)) * force[j]); - v_dot_f += velocity[j] * force[j]; - f_dot_f += force[j] * force[j]; + velocity[j] = velocity_new[j] + (h / (2 * mass)) * (force_prev[j] + force[j]); + + v_dot_f += velocity[j] * force[j]; + f_dot_f += force[j] * force[j]; force_prev[j] = force[j]; } if (v_dot_f <= 0) { - for(int j = 0; i < n_dofs_image; j++) { + for (int j = 0; i < n_dofs_image; j++) { velocity_new[j] = 0.0; } - } - else { - for(int j = 0; i < n_dofs_image; j++) { + } else { + for (int j = 0; i < n_dofs_image; j++) { velocity_new[j] = v_dot_f * force[j] / f_dot_f; } } - - } normalise_spins_C(y, n_images, n_dofs_image); diff --git a/fidimag/common/neb_method/nebm_integrators.h b/fidimag/common/neb_method/nebm_integrators.h index fac7e0f0..8f6da5b2 100644 --- a/fidimag/common/neb_method/nebm_integrators.h +++ b/fidimag/common/neb_method/nebm_integrators.h @@ -8,5 +8,4 @@ double step_Verlet_C(double *restrict forces, double mass, int n_images, int n_dofs_image, - double (* update_field) (double, double *) - ); + double (*update_field)(double, double *)); diff --git a/fidimag/common/neb_method/nebm_lib.c b/fidimag/common/neb_method/nebm_lib.c index 70d7ae5d..ec86b75c 100644 --- a/fidimag/common/neb_method/nebm_lib.c +++ b/fidimag/common/neb_method/nebm_lib.c @@ -2,7 +2,7 @@ #include "math.h" #include -void cross_product(double *restrict output, double * A, double * B){ +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 @@ -13,12 +13,12 @@ void cross_product(double *restrict output, double * A, double * B){ output[2] = A[0] * B[1] - A[1] * B[0]; } -double dot_product(double * A, double * B, int n){ +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++){ + for (int i = 0; i < n; i++) { dotp += A[i] * B[i]; } @@ -35,12 +35,14 @@ double compute_norm(double *restrict a, int n) { */ double norm = 0; - for(int i = 0; i < n; i++){ norm += a[i] * a[i]; } + for (int i = 0; i < n; i++) { + norm += a[i] * a[i]; + } norm = sqrt(norm); return norm; } -void normalise(double *restrict a, int n){ +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 @@ -55,56 +57,52 @@ void normalise(double *restrict a, int n){ length = compute_norm(a, n); - if (length > 0){ + if (length > 0) { length = 1.0 / length; } - for(int i = 0; i < n; i++){ + 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 *restrict y, int n_images, int n_dofs_image) { int i; // Index where the components of an image start in the *y array, int im_idx; - for(i = 1; i < n_images - 1; i++){ + for (i = 1; i < n_images - 1; i++) { im_idx = i * n_dofs_image; - double * y_i = &y[im_idx]; + double *y_i = &y[im_idx]; normalise(y_i, n_dofs_image); } } -void normalise_spins_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) { int i, j; // Index where the components of an image start in the *y array, int im_idx; - for(i = 1; i < n_images - 1; i++){ + for (i = 1; i < n_images - 1; i++) { im_idx = i * n_dofs_image; // Normalise per spin - for(j = 0; j < n_dofs_image; j+=3){ - double * y_i = &y[im_idx + j]; + for (j = 0; j < n_dofs_image; j += 3) { + double *y_i = &y[im_idx + j]; normalise(y_i, 3); } } } - /* ------------------------------------------------------------------------- */ - void compute_tangents_C(double *restrict tangents, double *restrict y, double *restrict energies, - int n_dofs_image, int n_images - ) { + int n_dofs_image, int n_images) { /* Compute the tangents for every degree of freedom of a NEBM band, * following the rules from Henkelman and Jonsson [Henkelman et al., @@ -167,31 +165,31 @@ 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_plus = malloc(n_dofs_image * sizeof(double)); t_minus = malloc(n_dofs_image * sizeof(double)); double deltaE_plus, deltaE_minus; double deltaE_MAX, deltaE_MIN; - for(i = 1; i < n_images - 1; i++){ + for (i = 1; i < n_images - 1; i++) { im_idx = i * (n_dofs_image); next_im_idx = (i + 1) * (n_dofs_image); prev_im_idx = (i - 1) * (n_dofs_image); // Tangents of the i-th image - double * t = &tangents[im_idx]; + double *t = &tangents[im_idx]; // Compute the t+ and t- vectors for the i-th image of the band, which // is given by the difference of the Y_i image with its neighbours - for(j = 0; j < n_dofs_image; j++){ - t_plus[j] = y[next_im_idx + j] - y[im_idx + j]; - t_minus[j] = y[im_idx + j] - y[prev_im_idx + j]; + for (j = 0; j < n_dofs_image; j++) { + t_plus[j] = y[next_im_idx + j] - y[im_idx + j]; + t_minus[j] = y[im_idx + j] - y[prev_im_idx + j]; } // Similarly, compute the energy differences - deltaE_plus = energies[i + 1] - energies[i]; - deltaE_minus = energies[i] - energies[i - 1]; + deltaE_plus = energies[i + 1] - energies[i]; + deltaE_minus = energies[i] - energies[i - 1]; /* Now we follow Henkelman and Jonsson rules [Henkelman et al., Journal * of Chemical Physics 113, 22 (2000)] for the tangent directions @@ -207,20 +205,22 @@ void compute_tangents_C(double *restrict tangents, double *restrict y, double *r /* ----------------------------------------------------------------- */ - if(deltaE_plus > 0 && deltaE_minus > 0) { - for(j = 0; j < n_dofs_image; j++) t[j] = t_plus[j]; + if (deltaE_plus > 0 && deltaE_minus > 0) { + for (j = 0; j < n_dofs_image; j++) + t[j] = t_plus[j]; } /* ----------------------------------------------------------------- */ - else if(deltaE_plus < 0 && deltaE_minus < 0) { - for(j = 0; j < n_dofs_image; j++) t[j] = t_minus[j]; + else if (deltaE_plus < 0 && deltaE_minus < 0) { + for (j = 0; j < n_dofs_image; j++) + t[j] = t_minus[j]; } /* ----------------------------------------------------------------- */ - else if((deltaE_plus < 0 && deltaE_minus > 0) || // A maximum - (deltaE_plus > 0 && deltaE_minus < 0)) { // A minimum + else if ((deltaE_plus < 0 && deltaE_minus > 0) || // A maximum + (deltaE_plus > 0 && deltaE_minus < 0)) { // A minimum /* According to the energy of the neighbouring images, the tangent * of the i-th image will be a combination of the differences wrt @@ -231,13 +231,13 @@ void compute_tangents_C(double *restrict tangents, double *restrict y, double *r deltaE_MIN = fmin(fabs(deltaE_plus), fabs(deltaE_minus)); if (energies[i + 1] > energies[i - 1]) { - for(j = 0; j < n_dofs_image; j++) { + for (j = 0; j < n_dofs_image; j++) { t[j] = t_plus[j] * deltaE_MAX + t_minus[j] * deltaE_MIN; } } else { - for(j = 0; j < n_dofs_image; j++) { + for (j = 0; j < n_dofs_image; j++) { t[j] = t_plus[j] * deltaE_MIN + t_minus[j] * deltaE_MAX; } } @@ -246,7 +246,8 @@ void compute_tangents_C(double *restrict tangents, double *restrict y, double *r // When energy has a constant slope else { - for(j = 0; j < n_dofs_image; j++) t[j] = t_plus[j] + t_minus[j]; + for (j = 0; j < n_dofs_image; j++) + t[j] = t_plus[j] + t_minus[j]; } /* ----------------------------------------------------------------- */ @@ -256,18 +257,16 @@ void compute_tangents_C(double *restrict tangents, double *restrict y, double *r free(t_minus); } // Close main function - /* ------------------------------------------------------------------------- */ 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 - ) { + double *restrict spring_force, + double *restrict y, + double *restrict tangents, + double *restrict k, + int n_images, + int n_dofs_image, + double *restrict distances) { /* Compute the spring force for every image of an energy band, which is * defined in the *y array. The *y array has (n_images * n_dofs_image) @@ -297,35 +296,32 @@ void compute_spring_force_C( double dY_plus_norm, dY_minus_norm; - for(i = 1; i < n_images - 1; i++){ + for (i = 1; i < n_images - 1; i++) { im_idx = i * (n_dofs_image); - double * sf = &spring_force[im_idx]; - double * t = &tangents[im_idx]; + double *sf = &spring_force[im_idx]; + double *t = &tangents[im_idx]; // Get the distances between the i-th image, Y_i, and its // neighbours, the (i+1)-th and (i-1)-th images. - dY_plus_norm = distances[i]; + dY_plus_norm = distances[i]; dY_minus_norm = distances[i - 1]; // Now compute the spring force - for(j = 0; j < n_dofs_image; j++) { + for (j = 0; j < n_dofs_image; j++) { sf[j] = k[i] * (dY_plus_norm - dY_minus_norm) * t[j]; } } - } - /* ------------------------------------------------------------------------- */ - void compute_effective_force_C(double *restrict G, double *restrict tangents, double *restrict gradientE, double *restrict spring_force, - int *restrict climbing_image, + int *restrict climbing_image, int n_images, int n_dofs_image) { @@ -357,30 +353,30 @@ void compute_effective_force_C(double *restrict G, int im_idx; double gradE_dot_t; - for(i = 1; i < n_images - 1; i++){ + for (i = 1; i < n_images - 1; i++) { im_idx = i * (n_dofs_image); - double * t = &tangents[im_idx]; - double * gradE = &gradientE[im_idx]; - double * sf = &spring_force[im_idx]; + double *t = &tangents[im_idx]; + double *gradE = &gradientE[im_idx]; + double *sf = &spring_force[im_idx]; gradE_dot_t = dot_product(gradE, t, n_dofs_image); - if(climbing_image[i] == 0) { - for(j = 0; j < n_dofs_image; j++) { + if (climbing_image[i] == 0) { + for (j = 0; j < n_dofs_image; j++) { G[im_idx + j] = -gradE[j] + gradE_dot_t * t[j] + sf[j]; } } // falling image with spring constant -> 0 - else if(climbing_image[i] == -1) { - for(j = 0; j < n_dofs_image; j++) { + else if (climbing_image[i] == -1) { + for (j = 0; j < n_dofs_image; j++) { G[im_idx + j] = -gradE[j]; } } // climbing image else { - for(j = 0; j < n_dofs_image; j++) { + for (j = 0; j < n_dofs_image; j++) { G[im_idx + j] = -gradE[j] + 2 * gradE_dot_t * t[j]; } } @@ -394,24 +390,23 @@ void compute_effective_force_C(double *restrict G, // degrees of freedom. This function updates the *distances and *path_distances // arrays. The *distances array saves the distances as: // [ |Y1 - Y0|, |Y2 - Y1|, ... ] -// Path distances are the total distances relative to the 0-th image +// 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, int n_images, int n_dofs_image, - double (* compute_distance)(double *, - double *, - int, - int *, - int), - int *restrict material, - int n_dofs_image_material - ) { + double (*compute_distance)(double *, + double *, + int, + int *, + int), + int *restrict material, + int n_dofs_image_material) { int i, im_idx, next_im_idx; path_distances[0] = 0.0; - for(i = 0; i < n_images - 1; i++){ + for (i = 0; i < n_images - 1; i++) { im_idx = i * (n_dofs_image); next_im_idx = (i + 1) * (n_dofs_image); @@ -419,8 +414,7 @@ void compute_image_distances(double *restrict distances, distances[i] = compute_distance(&y[next_im_idx], &y[im_idx], n_dofs_image, material, - n_dofs_image_material - ); + n_dofs_image_material); path_distances[i + 1] = path_distances[i] + distances[i]; } } @@ -428,8 +422,7 @@ void compute_image_distances(double *restrict distances, // ---------------------------------------------------------------------------- void project_vector_C(double *restrict vector, double *restrict y_i, - int n_dofs_image - ){ + int n_dofs_image) { /* Project a vector into the space of the y_i image in an energy band. We * assume that *image and *y_i have the same length. The vector and the @@ -459,20 +452,19 @@ void project_vector_C(double *restrict vector, double *restrict y_i, int j; double v_dot_m_i = 0; - for(j = 0; j < n_dofs_image; j++) { + for (j = 0; j < n_dofs_image; j++) { // Every 3 components of the *vector and *y_i array, compute the // 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_product(&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, - int n_images, int n_dofs_image - ){ + int n_images, int n_dofs_image) { /* Compute the projections of the vectors in the *vector array into the * space formed by the spin vectors in the array *y @@ -506,23 +498,21 @@ void project_images_C(double *restrict vector, double *restrict y, // Index where the components of an image start in the *y array, int im_idx; - for(i = 1; i < n_images - 1; i++){ + for (i = 1; i < n_images - 1; i++) { im_idx = i * (n_dofs_image); - double * v = &vector[im_idx]; - double * y_i = &y[im_idx]; + double *v = &vector[im_idx]; + double *y_i = &y[im_idx]; project_vector_C(v, y_i, 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 - ) { + int *restrict material, int n_dofs_image_material) { /* Compute the distance between two images, A and B, discarding the sites * without material @@ -535,9 +525,9 @@ double compute_distance_cartesian(double *restrict A, double *restrict B, int n_ double distance; int j = 0; - for(int i = 0; i < n_dofs_image; i++) { + for (int i = 0; i < n_dofs_image; i++) { if (material[i] > 0) { - A_minus_B[j] = A[i] - B[i]; + A_minus_B[j] = A[i] - B[i]; j += 1; } } @@ -570,19 +560,18 @@ 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, - int n_dofs_image - ) { + int n_dofs_image) { int n_spins = n_dofs_image / 3; - for(int i = 0; i < n_spins; i++){ - int j = 3 * i; + for (int i = 0; i < n_spins; i++) { + int j = 3 * i; - if (pins[i] > 0){ + if (pins[i] > 0) { dYdt[j] = 0; - dYdt[j + 1] = 0; - dYdt[j + 2] = 0; - continue; - } + dYdt[j + 1] = 0; + dYdt[j + 2] = 0; + 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]; @@ -592,11 +581,11 @@ void compute_dYdt(double *restrict Y, double *restrict G, double *restrict dYdt, dYdt[j + 2] = Y_dot_Y * G[j + 2] - Y_dot_G * Y[j + 2]; // Correction factor to rescale the spin length at every iteration step - double c = 6 * sqrt(dYdt[j] * dYdt[j] + + double c = 6 * sqrt(dYdt[j] * dYdt[j] + dYdt[j + 1] * dYdt[j + 1] + dYdt[j + 2] * dYdt[j + 2]); - dYdt[j] += c * (1 - Y_dot_Y) * Y[j]; + dYdt[j] += c * (1 - Y_dot_Y) * Y[j]; dYdt[j + 1] += c * (1 - Y_dot_Y) * Y[j + 1]; dYdt[j + 2] += c * (1 - Y_dot_Y) * Y[j + 2]; } @@ -604,9 +593,9 @@ 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, int n_images, int n_dofs_image) { - #pragma omp parallel for schedule(static) - for(int i = 1; i < n_images - 1; i++){ - //printf(""); +#pragma omp parallel for schedule(static) + for (int i = 1; i < n_images - 1; i++) { + // printf(""); int j = i * n_dofs_image; compute_dYdt(&y[j], &G[j], &dYdt[j], &pins[0], n_dofs_image); } @@ -619,19 +608,18 @@ void compute_dYdt_C(double *restrict y, double *restrict G, double *restrict dYd void compute_dYdt_nc(double *restrict Y, double *restrict G, double *restrict dYdt, int *restrict pins, - int n_dofs_image - ) { + int n_dofs_image) { int n_spins = n_dofs_image / 3; - for(int i = 0; i < n_spins; i++){ + for (int i = 0; i < n_spins; i++) { int j = 3 * i; - if (pins[i] > 0){ + if (pins[i] > 0) { dYdt[j] = 0; - dYdt[j + 1] = 0; - dYdt[j + 2] = 0; - continue; - } + dYdt[j + 1] = 0; + dYdt[j + 2] = 0; + 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]; @@ -645,9 +633,9 @@ void compute_dYdt_nc(double *restrict Y, double *restrict G, double *restrict dY void compute_dYdt_nc_C(double *restrict y, double *restrict G, double *restrict dYdt, int *restrict pins, int n_images, int n_dofs_image) { - #pragma omp parallel for schedule(static) - for(int i = 1; i < n_images - 1; i++){ - //printf(""); +#pragma omp parallel for schedule(static) + for (int i = 1; i < n_images - 1; i++) { + // printf(""); int j = i * n_dofs_image; compute_dYdt_nc(&y[j], &G[j], &dYdt[j], &pins[0], n_dofs_image); } diff --git a/fidimag/common/neb_method/nebm_lib.h b/fidimag/common/neb_method/nebm_lib.h index 6ebd9699..d666a373 100644 --- a/fidimag/common/neb_method/nebm_lib.h +++ b/fidimag/common/neb_method/nebm_lib.h @@ -1,67 +1,59 @@ #include "math.h" #define WIDE_PI 3.1415926535897932384626433832795L +void cross_product(double *output, double *A, double *B); -void cross_product(double * output, double * A, double * B); - -double dot_product(double * A, double * B, int n); +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 - ); + 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 - ); + 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 *restrict climbing_image, int n_images, - int n_dofs_image - ); + 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 - ); + 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 - ); + int n_images, int n_dofs_image); void project_vector_C(double *restrict vector, double *restrict y, - int n_dofs_image - ); + 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 - ); + int *restrict material, int n_dofs_image_material); void compute_dYdt_C(double *restrict y, double *restrict G, double *restrict dYdt, int *restrict pins, diff --git a/fidimag/common/neb_method/nebm_spherical_lib.c b/fidimag/common/neb_method/nebm_spherical_lib.c index 7ef9641f..7b588d74 100644 --- a/fidimag/common/neb_method/nebm_spherical_lib.c +++ b/fidimag/common/neb_method/nebm_spherical_lib.c @@ -1,25 +1,25 @@ #include "nebm_spherical_lib.h" -#include "nebm_lib.h" #include "math.h" +#include "nebm_lib.h" double compute_norm_spherical(double *restrict 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 @@ -28,26 +28,28 @@ double compute_norm_spherical(double *restrict a, int n, int scale) { */ double norm = 0; - double scale_db = (double) scale; + double scale_db = (double)scale; - for(int i = 0; i < n; i++){ + for (int i = 0; i < n; i++) { - if (a[i] > WIDE_PI){ + if (a[i] > WIDE_PI) { a[i] = 2 * WIDE_PI - a[i]; - } else if(a[i] < -WIDE_PI){ + } else if (a[i] < -WIDE_PI) { a[i] += 2 * WIDE_PI; } norm += a[i] * a[i]; } - if (scale == 0) norm = sqrt(norm); - else norm = sqrt(norm) / scale_db; + if (scale == 0) + norm = sqrt(norm); + else + norm = sqrt(norm) / scale_db; return norm; } -void normalise_spherical(double *restrict a, int n){ +void normalise_spherical(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 @@ -62,26 +64,25 @@ void normalise_spherical(double *restrict a, int n){ length = compute_norm_spherical(a, n, 0); - if (length > 0){ + if (length > 0) { length = 1.0 / length; } - for(int i = 0; i < n; i++){ + for (int i = 0; i < n; i++) { a[i] *= length; } } - -void normalise_images_spherical_C(double *restrict y, int n_images, int n_dofs_image){ +void normalise_images_spherical_C(double *restrict y, int n_images, int n_dofs_image) { int i; // Index where the components of an image start in the *y array, int im_idx; - for(i = 1; i < n_images - 1; i++){ + for (i = 1; i < n_images - 1; i++) { im_idx = i * n_dofs_image; - double * y_i = &y[im_idx]; + double *y_i = &y[im_idx]; normalise_spherical(y_i, n_dofs_image); } @@ -90,16 +91,15 @@ 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 - ) { + int *restrict material, int n_dofs_image_material) { double A_minus_B[n_dofs_image_material]; double distance; int j = 0; - for(int i = 0; i < n; i++) { + for (int i = 0; i < n; i++) { if (material[i] > 0) { - A_minus_B[j] = A[i] - B[i]; + A_minus_B[j] = A[i] - B[i]; j += 1; } } diff --git a/fidimag/common/neb_method/nebm_spherical_lib.h b/fidimag/common/neb_method/nebm_spherical_lib.h index 3e46ffce..db9d814a 100644 --- a/fidimag/common/neb_method/nebm_spherical_lib.h +++ b/fidimag/common/neb_method/nebm_spherical_lib.h @@ -6,6 +6,5 @@ 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 - ); +double compute_distance_spherical(double *restrict A, double *B, int n, + int *restrict material, int n_dofs_image_material); diff --git a/fidimag/micro/lib/anis.c b/fidimag/micro/lib/anis.c index a150042b..7f1b4ba9 100644 --- a/fidimag/micro/lib/anis.c +++ b/fidimag/micro/lib/anis.c @@ -1,17 +1,16 @@ #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) { -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; + int n = nx * ny * nz; - #pragma omp parallel for - for (int i = 0; i < n; i++) { - int j = 3 * i; +#pragma omp parallel for + for (int i = 0; i < n; i++) { + int j = 3 * i; - if (Ms_inv[i] == 0.0){ - field[j] = 0; + if (Ms_inv[i] == 0.0) { + field[j] = 0; field[j + 1] = 0; field[j + 2] = 0; energy[i] = 0; @@ -20,26 +19,23 @@ void compute_uniaxial_anis(double *restrict m, double *restrict field, double *r 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); - - } + 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 compute_uniaxial4_anis(double *restrict m, double *restrict field, double *restrict energy, double *restrict Ms_inv, + double *restrict K1, double *restrict K2, double *restrict axis, int nx, int ny, int nz) { -void compute_uniaxial4_anis(double *restrict m, double *restrict field, double *restrict energy, double *restrict Ms_inv, - double *restrict K1, double *restrict K2, double *restrict axis, int nx, int ny, int nz) { - int n = nx * ny * nz; // Follows calculation of OOMMF extension by Hans and Richard Boardman // http://www.soton.ac.uk/~fangohr/software/oxs_uniaxial4/download/uniaxialanisotropy4.cc - #pragma omp parallel for +#pragma omp parallel for for (int i = 0; i < n; i++) { int j = 3 * i; @@ -59,23 +55,22 @@ void compute_uniaxial4_anis(double *restrict m, double *restrict field, double * double m_dot_u = m[j] * axis[j] + m[j + 1] * axis[j + 1] + m[j + 2] * axis[j + 2]; if (k1 <= 0) { - field[j + 0] = (field_mult1*m_dot_u) * axis[j + 0] + (field_mult2 * m_dot_u*m_dot_u*m_dot_u) * axis[j + 0]; - field[j + 1] = (field_mult1*m_dot_u) * axis[j + 1] + (field_mult2 * m_dot_u*m_dot_u*m_dot_u) * axis[j + 1]; - field[j + 2] = (field_mult1*m_dot_u) * axis[j + 2] + (field_mult2 * m_dot_u*m_dot_u*m_dot_u) * axis[j + 2]; - energy[i] = -k1*m_dot_u*m_dot_u - k2*m_dot_u*m_dot_u*m_dot_u*m_dot_u; + field[j + 0] = (field_mult1 * m_dot_u) * axis[j + 0] + (field_mult2 * m_dot_u * m_dot_u * m_dot_u) * axis[j + 0]; + field[j + 1] = (field_mult1 * m_dot_u) * axis[j + 1] + (field_mult2 * m_dot_u * m_dot_u * m_dot_u) * axis[j + 1]; + field[j + 2] = (field_mult1 * m_dot_u) * axis[j + 2] + (field_mult2 * m_dot_u * m_dot_u * m_dot_u) * axis[j + 2]; + energy[i] = -k1 * m_dot_u * m_dot_u - k2 * m_dot_u * m_dot_u * m_dot_u * m_dot_u; } else { double u_x_m[3]; - u_x_m[0] = cross_x(axis[j], axis[j+1], axis[j+2], m[j], m[j+1], m[j+2]); - u_x_m[1] = cross_y(axis[j], axis[j+1], axis[j+2], m[j], m[j+1], m[j+2]); - u_x_m[2] = cross_z(axis[j], axis[j+1], axis[j+2], m[j], m[j+1], m[j+2]); - double u_x_m_mag2 = u_x_m[1]*u_x_m[1] + u_x_m[1]*u_x_m[1] + u_x_m[2]*u_x_m[2]; - field[j + 0] = (field_mult1*m_dot_u) * axis[j + 0] + (field_mult2*m_dot_u*m_dot_u*m_dot_u) * axis[j + 0]; - field[j + 1] = (field_mult1*m_dot_u) * axis[j + 1] + (field_mult2*m_dot_u*m_dot_u*m_dot_u) * axis[j + 1]; - field[j + 2] = (field_mult1*m_dot_u) * axis[j + 2] + (field_mult2*m_dot_u*m_dot_u*m_dot_u) * axis[j + 2]; - energy[i] = (k1 + 2*k2)*u_x_m_mag2 - k2*u_x_m_mag2*u_x_m_mag2; + u_x_m[0] = cross_x(axis[j], axis[j + 1], axis[j + 2], m[j], m[j + 1], m[j + 2]); + u_x_m[1] = cross_y(axis[j], axis[j + 1], axis[j + 2], m[j], m[j + 1], m[j + 2]); + u_x_m[2] = cross_z(axis[j], axis[j + 1], axis[j + 2], m[j], m[j + 1], m[j + 2]); + double u_x_m_mag2 = u_x_m[1] * u_x_m[1] + u_x_m[1] * u_x_m[1] + u_x_m[2] * u_x_m[2]; + field[j + 0] = (field_mult1 * m_dot_u) * axis[j + 0] + (field_mult2 * m_dot_u * m_dot_u * m_dot_u) * axis[j + 0]; + field[j + 1] = (field_mult1 * m_dot_u) * axis[j + 1] + (field_mult2 * m_dot_u * m_dot_u * m_dot_u) * axis[j + 1]; + field[j + 2] = (field_mult1 * m_dot_u) * axis[j + 2] + (field_mult2 * m_dot_u * m_dot_u * m_dot_u) * axis[j + 2]; + energy[i] = (k1 + 2 * k2) * u_x_m_mag2 - k2 * u_x_m_mag2 * u_x_m_mag2; } } - } \ No newline at end of file diff --git a/fidimag/micro/lib/dmi.c b/fidimag/micro/lib/dmi.c index bdb1e7b3..10ed6175 100644 --- a/fidimag/micro/lib/dmi.c +++ b/fidimag/micro/lib/dmi.c @@ -1,7 +1,7 @@ #include "micro_clib.h" -#include -#include #include +#include +#include /* Functions to Compute the micromagnetic Dzyaloshinskii Moriya interaction * field and energy using the @@ -151,23 +151,23 @@ * two DMIs: * * dmi_vector = [ Dx(-x) Dy(-x) Dz(-x) -- DMI vector 1 per every NN - * Dx(+x) Dy(+x) Dz(+x) - * Dx(-y) Dy(-y) Dz(-y) - * ... + * Dx(+x) Dy(+x) Dz(+x) + * Dx(-y) Dy(-y) Dz(-y) + * ... * Dx(+z) Dy(+z) Dz(+z) -- * Dx(-x) Dy(-x) Dz(-x) -- DMI vector 2 per every NN - * Dx(+x) Dy(+x) Dz(+x) - * Dx(-y) Dy(-y) Dz(-y) - * ... - * Dx(+z) Dy(+z) Dz(+z) + * Dx(+x) Dy(+x) Dz(+x) + * Dx(-y) Dy(-y) Dz(-y) + * ... + * Dx(+z) Dy(+z) Dz(+z) * ] * * and DMI constants are passed inerspersed in the D array per every node - * + * * D = [ D1_0 D2_0 D1_1 D2_1 ... D1_N D2_N] - * + * * where D1_i is the 1st DMI constant of the i-mesh node - * + * * ------------------------------------------------------------------------ * INTERFACIAL DMI * ------------------------------------------------------------------------ @@ -239,10 +239,8 @@ * continuum expression when doing the limit a_x a_y a_z --> 0 */ - // ---------------------------------------------------------------------------- - /* A function to compute any DMI given the corresponding DM vectors * obtained from discretising the Lifshitz invariants using the * finite differences: @@ -265,10 +263,10 @@ * the finite differences model. For DMIs defined in 2D * 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, + double *restrict D, int n_dmis, double *dmi_vector, double dx, double dy, double dz, int n, int *restrict ngbs) { @@ -280,10 +278,10 @@ void dmi_field(double *restrict m, double *restrict field, for (int i = 0; i < n; i++) { double DMIc; double fx = 0, fy = 0, fz = 0; - int idnm; // Index for the magnetisation matrix + 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){ + if (Ms_inv[i] == 0.0) { field[3 * i] = 0; field[3 * i + 1] = 0; field[3 * i + 2] = 0; @@ -302,7 +300,7 @@ void dmi_field(double *restrict m, double *restrict field, 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; @@ -312,19 +310,19 @@ void dmi_field(double *restrict m, double *restrict field, (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 - */ + 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 || + if (dmi_vector[ngbr_idx_D] != 0 || dmi_vector[ngbr_idx_D + 1] != 0 || - dmi_vector[ngbr_idx_D + 2] != 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) */ @@ -347,18 +345,16 @@ void dmi_field(double *restrict m, double *restrict field, dmi_vector[ngbr_idx_D + 1], dmi_vector[ngbr_idx_D + 2], m[idnm], m[idnm + 1], m[idnm + 2]); - } } - } // Close for loop through n of DMI constants - } // Close for loop through neighbours per mesh site + } // 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 * m[3 * i] + fy * m[3 * i + 1] - + fz * m[3 * i + 2]); + energy[i] = -0.5 * (fx * m[3 * i] + fy * m[3 * i + 1] + fz * m[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] = 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/fidimag/micro/lib/exch.c b/fidimag/micro/lib/exch.c index a12c4d69..c5d4c62b 100644 --- a/fidimag/micro/lib/exch.c +++ b/fidimag/micro/lib/exch.c @@ -1,8 +1,8 @@ #include "micro_clib.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) { + double *restrict Ms_inv, double A, double dx, double dy, double dz, + int n, int *restrict ngbs) { /* Compute the micromagnetic exchange field and energy using the * matrix of neighbouring spins and a second order approximation @@ -83,36 +83,36 @@ void compute_exch_field_micro(double *restrict m, double *restrict field, double */ /* Define the coefficients */ - double ax = 2 * A / (dx * dx); + double ax = 2 * A / (dx * dx); double ay = 2 * A / (dy * dy); double az = 2 * A / (dz * dz); /* Here we iterate through every mesh node */ - #pragma omp parallel for - for (int i = 0; i < n; i++) { - double fx = 0, fy = 0, fz = 0; - int idnm = 0; // Index for the magnetisation matrix - int idn = 6 * i; // index for the neighbours +#pragma omp parallel for + for (int i = 0; i < n; i++) { + 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; - } + 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) { + if (ngbs[idn + j] >= 0) { /* Magnetisation of the neighbouring spin since ngbs gives * the neighbour's index */ - idnm = 3 * ngbs[idn + j]; + idnm = 3 * ngbs[idn + j]; /* Check that the magnetisation of the neighbouring spin * is larger than zero */ - if (Ms_inv[ngbs[idn + j]] > 0){ + if (Ms_inv[ngbs[idn + j]] > 0) { /* Neighbours in the -x and +x directions * giving: ( m[i-x] - m[i] ) + ( m[i+x] - m[i] ) @@ -128,45 +128,44 @@ void compute_exch_field_micro(double *restrict m, double *restrict field, double * This same applies for the other directions */ if (j == 0 || j == 1) { - fx += ax * (m[idnm] - m[3 * i]); + fx += ax * (m[idnm] - m[3 * i]); fy += ax * (m[idnm + 1] - m[3 * i + 1]); fz += ax * (m[idnm + 2] - m[3 * i + 2]); } /* Neighbours in the -y and +y directions */ else if (j == 2 || j == 3) { - fx += ay * (m[idnm] - m[3 * i]); + fx += ay * (m[idnm] - m[3 * i]); fy += ay * (m[idnm + 1] - m[3 * i + 1]); fz += ay * (m[idnm + 2] - m[3 * i + 2]); } /* Neighbours in the -z and +z directions */ else if (j == 4 || j == 5) { - fx += az * (m[idnm] - m[3 * i]); + fx += az * (m[idnm] - m[3 * i]); fy += az * (m[idnm + 1] - m[3 * i + 1]); fz += az * (m[idnm + 2] - m[3 * i + 2]); + } else { + continue; } - else { - continue; } } } } /* Energy as: (-mu0 * Ms / 2) * [ H_ex * m ] */ - energy[i] = -0.5 * (fx * m[3 * i] + fy * m[3 * i + 1] - + fz * m[3 * i + 2]); + energy[i] = -0.5 * (fx * m[3 * i] + fy * m[3 * i + 1] + fz * m[3 * i + 2]); /* Update the field H_ex which has the same structure than *m */ - field[3 * i] = fx * Ms_inv[i] * MU0_INV; + 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; } } -inline int get_index(int nx, int ny, int i, int j, int k){ - return k * nx*ny + j * nx + i; +inline int get_index(int nx, int ny, int i, int j, int k) { + return k * nx * ny + j * nx + i; } 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){ + double sigma, int nx, double ny, double nz, int z_bottom, int z_top) { /* Compute the micromagnetic exchange field and energy using the * matrix of neighbouring spins and a second order approximation @@ -193,44 +192,42 @@ void compute_exch_field_rkky_micro(double *m, double *field, double *energy, dou * */ - int n = nx*ny*nz; - for (int i = 0; i < n; i++){ + int n = nx * ny * nz; + for (int i = 0; i < n; i++) { energy[i] = 0; - field[3*i]=0; - field[3*i+1]=0; - field[3*i+2]=0; - } - - #pragma omp parallel for - for (int i = 0; i < nx; i++) { - for (int j = 0; j < ny; j++){ - double mtx=0, mty=0, mtz=0; - double mbx=0, mby=0, mbz=0; - int id1 = get_index(nx,ny, i, j, z_bottom); - int id2 = get_index(nx,ny, i, j, z_top); - mtx = m[3*id2]; - mty = m[3*id2+1]; - mtz = m[3*id2+2]; - - mbx = m[3*id1]; - mby = m[3*id1+1]; - mbz = m[3*id1+2]; - - if (Ms_inv[id1] != 0.0){ - energy[id1] = sigma*(1-mtx*mbx-mty*mby-mtz*mbz); - field[3*id1] = sigma * mtx * Ms_inv[id1] * MU0_INV; - field[3*id1+1] = sigma * mty * Ms_inv[id1] * MU0_INV; - field[3*id1+2] = sigma * mtz * Ms_inv[id1] * MU0_INV; - } - - if (Ms_inv[id2] != 0.0){ - energy[id2] = sigma*(1-mtx*mbx-mty*mby-mtz*mbz); - field[3*id2] = sigma * mbx * Ms_inv[id2] * MU0_INV; - field[3*id2+1] = sigma * mby * Ms_inv[id2] * MU0_INV; - field[3*id2+2] = sigma * mbz * Ms_inv[id2] * MU0_INV; - } - } - } + field[3 * i] = 0; + field[3 * i + 1] = 0; + field[3 * i + 2] = 0; + } +#pragma omp parallel for + for (int i = 0; i < nx; i++) { + for (int j = 0; j < ny; j++) { + double mtx = 0, mty = 0, mtz = 0; + double mbx = 0, mby = 0, mbz = 0; + int id1 = get_index(nx, ny, i, j, z_bottom); + int id2 = get_index(nx, ny, i, j, z_top); + mtx = m[3 * id2]; + mty = m[3 * id2 + 1]; + mtz = m[3 * id2 + 2]; + + mbx = m[3 * id1]; + mby = m[3 * id1 + 1]; + mbz = m[3 * id1 + 2]; + + if (Ms_inv[id1] != 0.0) { + energy[id1] = sigma * (1 - mtx * mbx - mty * mby - mtz * mbz); + field[3 * id1] = sigma * mtx * Ms_inv[id1] * MU0_INV; + field[3 * id1 + 1] = sigma * mty * Ms_inv[id1] * MU0_INV; + field[3 * id1 + 2] = sigma * mtz * Ms_inv[id1] * MU0_INV; + } + if (Ms_inv[id2] != 0.0) { + energy[id2] = sigma * (1 - mtx * mbx - mty * mby - mtz * mbz); + field[3 * id2] = sigma * mbx * Ms_inv[id2] * MU0_INV; + field[3 * id2 + 1] = sigma * mby * Ms_inv[id2] * MU0_INV; + field[3 * id2 + 2] = sigma * mbz * Ms_inv[id2] * MU0_INV; + } + } + } } diff --git a/fidimag/micro/lib/micro_clib.h b/fidimag/micro/lib/micro_clib.h index 9a1bd2cf..f38bdf21 100644 --- a/fidimag/micro/lib/micro_clib.h +++ b/fidimag/micro/lib/micro_clib.h @@ -1,20 +1,20 @@ -#include +#include -#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; } + 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; } + 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; } + 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); + 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, @@ -23,13 +23,13 @@ void dmi_field(double *restrict m, double *restrict field, 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); + 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 *restrict Ku, double *restrict axis, int nx, int ny, int nz); -void compute_uniaxial4_anis(double *restrict m, double *restrict field, double *restrict energy, double *restrict Ms_inv, - double *restrict K1, double *restrict K2, double *restrict axis, int nx, int ny, int nz); +void compute_uniaxial4_anis(double *restrict m, double *restrict field, double *restrict energy, double *restrict Ms_inv, + double *restrict K1, double *restrict K2, 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/util.c b/fidimag/micro/lib/util.c index 6982fe3e..415a2569 100644 --- a/fidimag/micro/lib/util.c +++ b/fidimag/micro/lib/util.c @@ -2,10 +2,10 @@ // Compute: S \cdot (S_i \times S_j) inline double volume(double S[3], double Si[3], double Sj[3]) { - double tx = S[0] * (-Si[2] * Sj[1] + Si[1] * Sj[2]); - double ty = S[1] * (Si[2] * Sj[0] - Si[0] * Sj[2]); - double tz = S[2] * (-Si[1] * Sj[0] + Si[0] * Sj[1]); - return tx + ty + tz; + double tx = S[0] * (-Si[2] * Sj[1] + Si[1] * Sj[2]); + double ty = S[1] * (Si[2] * Sj[0] - Si[0] * Sj[2]); + double tz = S[2] * (-Si[1] * Sj[0] + Si[0] * Sj[1]); + return tx + ty + tz; } double skyrmion_number(double *restrict spin, double *restrict charge, @@ -35,7 +35,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 ... ] @@ -52,28 +52,28 @@ double skyrmion_number(double *restrict spin, double *restrict charge, * in the -y direction, for example */ - int i, j; - int index, id; + int i, j; + int index, id; - double sum = 0; + double sum = 0; - double S[3]; + double S[3]; - int nxy = nx * ny; + int nxy = nx * ny; /* 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 - } - for (i = 0; i < nxy; i++) { + for (i = 0; i < 12; i++) { + S_nn[i] = 0; // we have to set S_nn to zeros manually + } + + 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; @@ -84,22 +84,18 @@ double skyrmion_number(double *restrict spin, double *restrict charge, for (j = 0; j < 4; j++) { if (ngbs[id_nn + j] > 0) { id = 3 * ngbs[id_nn + j]; - S_nn[3 * j ] = spin[id]; + S_nn[3 * j] = spin[id]; S_nn[3 * j + 1] = spin[id + 1]; S_nn[3 * j + 2] = spin[id + 2]; } } - charge[i] = volume(S, &S_nn[3], &S_nn[9]) - + volume(S, &S_nn[0], &S_nn[6]) - - volume(S, &S_nn[0], &S_nn[9]) - - volume(S, &S_nn[3], &S_nn[6]); + charge[i] = volume(S, &S_nn[3], &S_nn[9]) + volume(S, &S_nn[0], &S_nn[6]) - volume(S, &S_nn[0], &S_nn[9]) - volume(S, &S_nn[3], &S_nn[6]); charge[i] /= (16 * WIDE_PI); sum += charge[i]; } - return sum; - + return sum; } diff --git a/tests/integrators_test.py b/tests/integrators_test.py index e95050f0..ff90dac2 100644 --- a/tests/integrators_test.py +++ b/tests/integrators_test.py @@ -33,7 +33,10 @@ def test_step(integrator, stepsize, debug=False): ys[i+1] = yp if not debug: assert 85 < ys[-1] < 100 - return ts, ys + + # Tests should return None (see PytestReturnNotNoneWarning) + # TODO: Compare solutions with expected one + # return ts, ys @pytest.mark.parametrize("integrator,stepsize_reported,stepsize_internal", [ ("euler", 0.2, 0.2), @@ -54,9 +57,11 @@ def test_step_integrator(integrator, stepsize_reported, stepsize_internal): ts[i+1] = integrator.t ys[i+1] = integrator.y assert 85 < ys[-1] < 100 - return ts, ys + + # return ts, ys +@pytest.mark.skip(reason="Not using assert") def test_scipy_integrator(): y_true = lambda t: np.sin(t) + t f = lambda t, y: np.cos(t) + 1 # derivative of f we'll use for integration @@ -73,13 +78,14 @@ def test_scipy_integrator(): for i, t in enumerate(ts[1:]): od.run_until(t) ts[i+1] = od.t - ys[i+1] = od.y + ys[i+1] = od.y[0] print("t = {}".format(t)) print(len(od.internal_timesteps)) print(od.rhs_evals) - return y_true, ts, ys, od.internal_timesteps + + # return y_true, ts, ys, od.internal_timesteps if __name__ == "__main__":