Commit final thesis version of FMM/BH in Fidimag

Author: rpep

fidimag/atomistic/demag.py

@@ -103,7 +103,13 @@ def compute_energy(self):
 
 
 class DemagFMM(Energy):
-    def __init__(self, order, ncrit, theta, name="DemagFMM"):
+    def __init__(self, order, ncrit, theta, name="DemagFMM", type='fmm'):
+        self.type = type
+        if self.type == 'fmm':
+            self._type = 0
+        elif self.type == 'bh':
+            self._type = 1
+
         self.name = name
         assert order > 0, "Order must be 1 or higher"
         assert order < 11, "Order bust be < 11"
@@ -121,11 +127,12 @@ def setup(self, mesh, spin, mu_s, mu_s_inv):
                            self.order,
                            self.coords,
                            self.m_temp,
-                           self.mu_s,
+                           self.mu_s, self._type
                            )
 
     def compute_field(self, t=0, spin=None):
+        print('computing field')
         self.m_temp[:] = spin if spin is not None else self.spin
         self.fmm.compute_field(self.field)
-        self.field *= 1e-7
+        self.field *= -1e-7
         return self.field

fidimag/atomistic/fmmlib/calculate.cpp

@@ -6,60 +6,44 @@
 #include <stack>
 #include <cmath>
 
-void P2P(double x, double y, double z, double mux, double muy, double muz, double *F) {
-  double R2 = x*x + y*y + z*z;
-  double R3 = R2*sqrt(R2);
-  double R5 = R3*R2;
-  double mu_dot_r = mux*x + muy*y + muz*z;
-  #pragma omp atomic
-  F[0] += (3*mu_dot_r * x / R5 - mux / R3);
-  #pragma omp atomic
-  F[1] += (3*mu_dot_r * y / R5 - muy / R3);
-  #pragma omp atomic
-  F[2] += (3*mu_dot_r * z / R5 - muz / R3);
-}
-
-void P2P_noatomic(double x, double y, double z, double mux, double muy, double muz, double *F) {
-  double R2 = x*x + y*y + z*z;
-  double R3 = R2*sqrt(R2);
-  double R5 = R3*R2;
-  double mu_dot_r = mux*x + muy*y + muz*z;
-  F[0] += (3*mu_dot_r * x / R5 - mux / R3);
-  F[1] += (3*mu_dot_r * y / R5 - muy / R3);
-  F[2] += (3*mu_dot_r * z / R5 - muz / R3);
+void M_sanity_check(const std::vector<Cell> &cells) {
+	double M0 = 0;
+	for(size_t c = 1; c < cells.size(); c++) {
+      if (cells[c].nchild == 0) {
+		    M0 += cells[c].M[0];
+	    }
+  }
+	std::cout << "Cell 0 has M0 = " << cells[0].M[0] << std::endl;
+	std::cout << "Other cells   = " << M0 << std::endl;
+  if (std::abs((cells[0].M[0] - M0)/M0) > 10e-10) {
+    throw std::runtime_error("M0 sanity check failed");
+  }
 }
 
 
-void P2P_Cells(size_t A, size_t B, std::vector<Cell> &cells,
-  std::vector<Particle> &particles, double *F) {
+void P2P_Cells(size_t A, size_t B, std::vector<Cell> &cells, std::vector<Particle> &particles, double *F) {
   // A - target
   // B - source
   for (size_t p1 = 0; p1 < cells[A].nleaf; p1++) {
-    double F_p[3] = {0.0};
+    //double F_p[FMMGEN_OUTPUTSIZE] = {0.0};
     size_t l1 = cells[A].leaf[p1];
     for (size_t p2 = 0; p2 < cells[B].nleaf; p2++) {
       size_t l2 = cells[B].leaf[p2];
       if (l2 != l1) {
-      	double dx = particles[l1].r[0] - particles[l2].r[0];
-      	double dy = particles[l1].r[1] - particles[l2].r[1];
-      	double dz = particles[l1].r[2] - particles[l2].r[2];
-      	P2P(dx, dy, dz, particles[l2].mu[0], particles[l2].mu[1], particles[l2].mu[2], F_p);
+      	double dx = (particles[l1].r[0] - particles[l2].r[0]);
+      	double dy = (particles[l1].r[1] - particles[l2].r[1]);
+      	double dz = (particles[l1].r[2] - particles[l2].r[2]);
+      	P2P(dx, dy, dz, particles[l2].S, &F[FMMGEN_OUTPUTSIZE*l1]);
       }
     }
-    #pragma omp atomic
-    F[3 * l1 + 0] += F_p[0];
-    #pragma omp atomic
-    F[3 * l1 + 1] += F_p[1];
-    #pragma omp atomic
-    F[3 * l1 + 2] += F_p[2];
   }
 }
 
 void interact_dehnen(size_t A, size_t B, std::vector<Cell> &cells, std::vector<Particle> &particles,
 		     double theta, size_t order, size_t ncrit, double *F) {
-  double dx = cells[A].x - cells[B].x;
-  double dy = cells[A].y - cells[B].y;
-  double dz = cells[A].z - cells[B].z;
+  double dx = (cells[A].x - cells[B].x);
+  double dy = (cells[A].y - cells[B].y);
+  double dz = (cells[A].z - cells[B].z);
   double R = sqrt(dx*dx + dy*dy + dz*dz);
 
   if (R*theta > (cells[A].rmax + cells[B].rmax)) {
@@ -97,7 +81,7 @@ void interact_dehnen(size_t A, size_t B, std::vector<Cell> &cells, std::vector<P
 void interact_dehnen_lazy(const size_t A, const size_t B,
                           const std::vector<Cell> &cells,
                           const std::vector<Particle> &particles,
-			              const double theta, const size_t order,
+			                    const double theta, const size_t order,
                           const size_t ncrit,
                           std::vector<std::pair<size_t, size_t>> &M2L_list,
                           std::vector<std::pair<size_t, size_t>> &P2P_list) {
@@ -150,32 +134,30 @@ void interact_dehnen_lazy(const size_t A, const size_t B,
 
 void evaluate_P2M(std::vector<Particle> &particles, std::vector<Cell> &cells,
               size_t cell, size_t ncrit, size_t exporder) {
-  // if (cells[cell].nleaf >= ncrit) {
-  //   for (size_t octant = 0; octant < 8; octant++) {
-  //     if (cells[cell].nchild & (1 << octant)) {
-  //     evaluate_P2M(particles, cells, cells[cell].child[octant], ncrit, exporder);
-  //     }
-  //   }
-  // }
-  // else {
-  double *M = new double[Nterms(exporder+1)]();
   #pragma omp for
   for(size_t c = 0; c < cells.size(); c++) {
+    //std::cout << "Cell " << c << std::endl;
+    //std::cout << "  Msize = " << Msize(exporder, FMMGEN_SOURCEORDER) << std::endl;
+    size_t msize = Msize(exporder, FMMGEN_SOURCEORDER);
+    double *M = new double[msize]();
+
     if (cells[c].nleaf < ncrit) {
-    for(size_t i = 0; i < cells[c].nleaf; i++) {
-      size_t l = cells[c].leaf[i];
-      M[0] = particles[l].mu[0];
-      M[1] = particles[l].mu[1];
-      M[2] = particles[l].mu[2];
-      // std::cout << "mu[" << l << "] = " << particles[l].mu[0] << std::endl;
-      double dx = (particles[l].r[0] - cells[c].x);
-      double dy = (particles[l].r[1] - cells[c].y);
-      double dz = (particles[l].r[2] - cells[c].z);
-      M2M(-dx, -dy, -dz, M, cells[c].M, exporder);
-    }
+      for(size_t i = 0; i < cells[c].nleaf; i++) {
+        size_t l = cells[c].leaf[i];
+        // Walter dehnen's definition:
+        // (-1)^m / m! (x_a - z_a)^m
+        double dx = (cells[c].x - particles[l].r[0]);
+        double dy = (cells[c].y - particles[l].r[1]);
+        double dz = (cells[c].z - particles[l].r[2]);
+        for(int k = 0; k < FMMGEN_SOURCESIZE; k++) {
+          // std::cout << particles[l].S[k] << std::endl;
+          M[k] = particles[l].S[k];
+        }
+        M2M(dx, dy, dz, M, cells[c].M, exporder);
+      }
    }
+   delete[] M;
   }
-  delete[] M;
 }
 
 void evaluate_M2M(std::vector<Particle> &particles, std::vector<Cell> &cells,
@@ -189,12 +171,17 @@ void evaluate_M2M(std::vector<Particle> &particles, std::vector<Cell> &cells,
   by iterating backwards through the nodes because
   of the way the tree is constructed.
   */
-  #pragma omp for
-  for (size_t c = 1; c < cells.size(); c++) {
+  // #pragma omp for schedule(dynamic)
+  // Can't currently go up the tree in parallel.
+  // Needs to be recursive or summing not correct.
+
+  // Dehnen definition:
+  // M_m(z_p) = (z_p - z_c)^n / n! M_{m - n}
+  for (size_t c = cells.size() - 1; c > 0; c--) {
     size_t p = cells[c].parent;
-    double dx = cells[p].x - cells[c].x;
-    double dy = cells[p].y - cells[c].y;
-    double dz = cells[p].z - cells[c].z;
+    double dx = (cells[p].x - cells[c].x);
+    double dy = (cells[p].y - cells[c].y);
+    double dz = (cells[p].z - cells[c].z);
     M2M(dx, dy, dz, cells[c].M, cells[p].M, exporder);
   }
 }
@@ -206,6 +193,11 @@ void evaluate_M2L_lazy(std::vector<Cell> &cells,
     for(size_t i = 0; i < M2L_list.size(); i++) {
     	size_t B = M2L_list[i].first;
     	size_t A = M2L_list[i].second;
+      // Dehnen definition:
+      // F_n(z_B) = M_m(z_A) * D_{n+m} (z_B - z_A)
+
+      // So here, we've reversed the order:
+      // F_n(z_A) = M_m(z_B) * D_{n+m} (z_A - z_B)
     	double dx = cells[A].x - cells[B].x;
     	double dy = cells[A].y - cells[B].y;
     	double dz = cells[A].z - cells[B].z;
@@ -225,16 +217,17 @@ void evaluate_P2P_lazy(std::vector<Cell> &cells, std::vector<Particle> &particle
 
 
 void evaluate_L2L(std::vector<Cell> &cells, size_t exporder) {
-  #pragma omp for // schedule(runtime)
-  for (size_t i = 0; i < cells.size(); i++) {
+  // Can't currently go down the tree in parallel!
+  // needs to be recursive or summing not correct.
+  for (size_t p = 0; p < cells.size(); p++) {
     for (int octant = 0; octant < 8; octant++) {
-      if (cells[i].nchild & (1 << octant)) {
-        // for child in cell i
-        size_t c = cells[i].child[octant];
-        double dx = cells[c].x - cells[i].x;
-        double dy = cells[c].y - cells[i].y;
-        double dz = cells[c].z - cells[i].z;
-        L2L(dx, dy, dz, cells[i].L, cells[c].L, exporder);
+      if (cells[p].nchild & (1 << octant)) {
+        // for child c in cell p
+        size_t c = cells[p].child[octant];
+        double dx = cells[c].x - cells[p].x;
+        double dy = cells[c].y - cells[p].y;
+        double dz = cells[c].z - cells[p].z;
+        L2L(dx, dy, dz, cells[p].L, cells[c].L, exporder);
       }
     }
   }
@@ -245,36 +238,77 @@ void evaluate_L2P(std::vector<Particle> &particles, std::vector<Cell> &cells,
   #pragma omp for schedule(runtime)
   for (size_t i = 0; i < cells.size(); i++) {
     if (cells[i].nleaf < ncrit) {
-      // std::cout << "cell " << i << " is a leaf " << std::endl;
       for (size_t p = 0; p < cells[i].nleaf; p++) {
 	    size_t k = cells[i].leaf[p];
-	    // std::cout << "L2P from " << i << " to particle " << k << std::endl;
         double dx = particles[k].r[0] - cells[i].x;
         double dy = particles[k].r[1] - cells[i].y;
         double dz = particles[k].r[2] - cells[i].z;
-	    double Fv[3] = {0.0};
-        L2P(dx, dy, dz, cells[i].L, Fv, exporder);
-    	F[3*k+0] -= Fv[0];
-    	F[3*k+1] -= Fv[1];
-    	F[3*k+2] -= Fv[2];
+        L2P(dx, dy, dz, cells[i].L, &F[FMMGEN_OUTPUTSIZE*k], exporder);
       }
     }
   }
 }
 
 void evaluate_direct(std::vector<Particle> &particles, double *F, size_t n) {
- #pragma omp parallel for schedule(runtime)
+  #pragma omp parallel for schedule(runtime)
   for (size_t i = 0; i < n; i++) {
     for (size_t j = 0; j < n; j++) {
-        if (i != j) {
-        	double dx = particles[i].r[0] - particles[j].r[0];
-        	double dy = particles[i].r[1] - particles[j].r[1];
-        	double dz = particles[i].r[2] - particles[j].r[2];
-        	double mux = particles[j].mu[0];
-        	double muy = particles[j].mu[1];
-        	double muz = particles[j].mu[2];
-        	P2P(dx, dy, dz, mux, muy, muz, &F[3*i]);
+      if (i != j) {
+      	double dx = particles[i].r[0] - particles[j].r[0];
+      	double dy = particles[i].r[1] - particles[j].r[1];
+      	double dz = particles[i].r[2] - particles[j].r[2];
+      	P2P(dx, dy, dz, particles[j].S, &F[FMMGEN_OUTPUTSIZE*i]);
+      }
+    }
+  }
+}
+
+void evaluate_M2P_and_P2P(std::vector<Particle> &particles, unsigned int p, unsigned int i,
+  std::vector<Cell> &cells, double *F, unsigned int n_crit, double theta,
+  unsigned int exporder) {
+  // For particle i, start at cell p
+  double dx, dy, dz, r;
+  int c;
+  unsigned int j;
+  // If cell p is not a leaf cell:
+  if (cells[p].nleaf >= n_crit) {
+    // Iterate through it's children
+    for (unsigned int octant = 0; octant < 8; octant++) {
+      // If a child exists in a given octant:
+      if (cells[p].nchild & (1 << octant)) {
+        // Get the child's index
+        c = cells[p].child[octant];
+        // Calculate the distance from the particle to the child cell
+        dx = particles[i].r[0] - cells[c].x;
+        dy = particles[i].r[1] - cells[c].y;
+        dz = particles[i].r[2] - cells[c].z;
+        r = sqrt(dx*dx + dy*dy + dz*dz);
+        // Apply the Barnes-Hut criterion:
+        if (cells[c].r > theta * r) {
+            // If the cell is 'near':
+            evaluate_M2P_and_P2P(particles, c, i, cells, F, n_crit, theta, exporder);
+        }
+        else {
+            // If the cell is 'far', calculate the potential and field
+            // on the particle from the multipole expansion:
+            M2P(dx, dy, dz, cells[c].M, &F[FMMGEN_OUTPUTSIZE*i], exporder);
         }
       }
+    }
+  }
+  else {
+    // loop in leaf cell's particles
+    for(unsigned int l = 0; l < (cells[p].nleaf); l++) {
+      // Get the particle index:
+      j = cells[p].leaf[l];
+      if (i != j) {
+        // Calculate the interparticle distance:
+        dx = particles[i].r[0] - particles[j].r[0];
+        dy = particles[i].r[1] - particles[j].r[1];
+        dz = particles[i].r[2] - particles[j].r[2];
+        // Compute the field:
+        P2P(dx, dy, dz, particles[j].S, &F[FMMGEN_OUTPUTSIZE*i]);
+      }
+    }
   }
 }

fidimag/atomistic/fmmlib/calculate.hpp

@@ -14,7 +14,7 @@
 #include <iostream>
 #include<omp.h>
 
-void P2P(double x, double y, double z, double mux, double muy, double muz, double *F);
+void M_sanity_check(const std::vector<Cell> &cells);
 
 void evaluate_P2M(std::vector<Particle> &particles, std::vector<Cell> &cells,
 		  size_t cell, size_t ncrit, size_t exporder);
@@ -52,3 +52,7 @@ void evaluate_M2L_lazy(std::vector<Cell> &cells,
 
 void evaluate_P2P_lazy(std::vector<Cell> &cells, std::vector<Particle> &particles,
                     std::vector<std::pair<size_t, size_t>> &P2P_list, double *F);
+
+void evaluate_M2P_and_P2P(std::vector<Particle> &particles, unsigned int p, unsigned int i,
+   std::vector<Cell> &cells, double *F, unsigned int n_crit, double theta,
+   unsigned int exporder);

fidimag/atomistic/fmmlib/example.py

@@ -1,14 +1,15 @@
 import fmmgen
 
-order = 11
+order = 8
 source_order = 1
 cse = True
 atomic = True
 precision='double'
 fmmgen.generate_code(order, "operators",
                      precision=precision,
                      CSE=cse,
-                     cython_wrapper=False,
+                     cython=False,
+                     harmonic_derivs=True,
                      potential=False,
                      field=True,
                      source_order=source_order,