Skip to content

Commit

Permalink
Clang formatted C extensions. Updated integrators test
Browse files Browse the repository at this point in the history
  • Loading branch information
@davidcortesortuno
davidcortesortuno committed Feb 2, 2024
1 parent eeb9eba commit d0ac9fb
Show file tree
Hide file tree
Showing 23 changed files with 1,407 additions and 1,543 deletions.
855 changes: 414 additions & 441 deletions fidimag/common/dipolar/demag.c
View file

Large diffs are not rendered by default.

526 changes: 257 additions & 269 deletions fidimag/common/dipolar/demag_oommf.c
View file

Large diffs are not rendered by default.

104 changes: 46 additions & 58 deletions fidimag/common/dipolar/demagcoef.h
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -39,10 +39,6 @@ to copyright protection within the United States.
*/





/* FILE: demagcoef.h -*-Mode: c++-*-
*
* Demag coefficients.
Expand Down Expand Up @@ -73,106 +69,98 @@ to copyright protection within the United States.
*
*/

#include <math.h>

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(x<y) return 1;
if(x>y) 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
// from extra precision. Also, don't expect the compiler to respect
// order with respect to parentheses at high levels of optimization,
// i.e., write "u=x; u-=(y-corr)" as opposed to "u=x-(y-corr)".

double sum,corr,y,u,t,v,z,x,tmp;

sum=arr[0]; corr=0;
for(int i=1;i<n;i++) {
x=arr[i];
y=corr+x;
tmp=y-corr;
u=x-tmp;
t=y+sum;
tmp=t-sum;
v=y-tmp;
z=u+v;
sum=t+z;
tmp=sum-t;
corr=z-tmp;
double sum, corr, y, u, t, v, z, x, tmp;

sum = arr[0];
corr = 0;
for (int i = 1; i < n; i++) {
x = arr[i];
y = corr + x;
tmp = y - corr;
u = x - tmp;
t = y + sum;
tmp = t - sum;
v = y - tmp;
z = u + v;
sum = t + z;
tmp = sum - t;
corr = z - tmp;
}
return sum;
}

155 changes: 78 additions & 77 deletions fidimag/common/dipolar/dipolar.h
View file
Original file line number Diff line number Diff line change
@@ -1,77 +1,81 @@
#include<math.h>
#include<complex.h>
#include<fftw3.h>
//#include<omp.h>
#include <complex.h>
#include <fftw3.h>
#include <math.h>
// #include<omp.h>

#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);

Expand All @@ -82,36 +86,33 @@ 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;

void init_solver(ode_solver *s, double k_B, double theta, int nxyz, double dt, double gamma);
ode_solver *create_ode_plan(void);
void finalize_ode_plan(ode_solver *plan);
void run_step1(ode_solver *s, double *m, double *h, double *m_pred, double *T,
double *alpha, double *mu_s_inv, int *pins);
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);
Loading

0 comments on commit d0ac9fb

Please sign in to comment.