Skip to content
Snippets Groups Projects
Select Git revision
  • MPI_qsched
  • resource_reuse
  • task_timers
  • openMP_locks
  • fortran_bindings
  • fortran_with_timers
  • master default protected
  • paper_updates
  • aidan
  • ompss_fixes
  • paper_fixes
  • FMM
  • autotools-update
  • paper_simulated_memcpy
  • gpu_friendly_bh
  • unbiased_gravity
  • test_bh_2
  • cleanup
  • merge_from_svn
19 results
Blame Permalink
multipoles.h 9.96 KiB 


struct multipole {

#ifdef QUADRUPOLES
  /* Quadrupole */
  float M_200;
  float M_020;
  float M_002;
  float M_110;
  float M_101;
  float M_011;
#endif
  /* CoM */
  double M_100;
  double M_010;
  double M_001;
 
  /* Mass */
  float M_000;
};


struct fieldTensor {

#ifdef QUADRUPOLES
  float F_200;
  float F_020;
  float F_002;
  float F_110;
  float F_101;
  float F_011;
#endif
  float F_100;
  float F_010;
  float F_001;

  float F_000;
};


static inline void initFieldTensor(struct fieldTensor f) {

  f.F_000 = 0.f;

  f.F_100 = f.F_010 = f.F_001 = 0.f;

#ifdef QUADRUPOLES
  f.F_200 = f.F_020 = f.F_002 = 0.f;
  f.F_110 = f.F_101 = f.F_011 = 0.f;
#endif

}


/******************************
 * Derivatives of the potential
 * 
 * All assume inv_r = 1/sqrt(r_x^2 + r_y^2 + r_z^2)
 ******************************/

/* 0-th order */
static inline double D_000(double r_x, double r_y, double r_z, double inv_r) {  // phi(r)

   return  inv_r; 
}

/* 1-st order */
static inline double D_100(double r_x, double r_y, double r_z, double inv_r) {  // d/dx phi(r)

  return -r_x * inv_r * inv_r * inv_r; 
}  

static inline double D_010(double r_x, double r_y, double r_z, double inv_r) {  // d/dy phi(r)

  return -r_y * inv_r * inv_r * inv_r; 
}  

static inline double D_001(double r_x, double r_y, double r_z, double inv_r) {  // d/dz phi(r)

  return -r_z * inv_r * inv_r * inv_r; 
}

/* 2-nd order */
static inline double D_200(double r_x, double r_y, double r_z, double inv_r) {  // d2/dxdx phi(r)
  
  double inv_r2 = inv_r * inv_r;
  double inv_r3 = inv_r * inv_r2;
  double inv_r5 = inv_r3 * inv_r2;
  return 3. * r_x * r_x * inv_r5 - inv_r3;
}

static inline double D_020(double r_x, double r_y, double r_z, double inv_r) {  // d2/dydy phi(r)
  
  double inv_r2 = inv_r * inv_r;
  double inv_r3 = inv_r * inv_r2;
  double inv_r5 = inv_r3 * inv_r2;
  return 3. * r_y * r_y * inv_r5 - inv_r3;
}

static inline double D_002(double r_x, double r_y, double r_z, double inv_r) {  // d2/dzdz phi(r)
  
  double inv_r2 = inv_r * inv_r;
  double inv_r3 = inv_r * inv_r2;
  double inv_r5 = inv_r3 * inv_r2;
  return 3. * r_z * r_z * inv_r5 - inv_r3;
}

static inline double D_110(double r_x, double r_y, double r_z, double inv_r) {  // d2/dxdy phi(r)

  double inv_r2 = inv_r * inv_r;
  double inv_r5 = inv_r2 * inv_r2 * inv_r;
  return 3. * r_x * r_y * inv_r5;
}

static inline double D_101(double r_x, double r_y, double r_z, double inv_r) {  // d2/dxdz phi(r)

  double inv_r2 = inv_r * inv_r;
  double inv_r5 = inv_r2 * inv_r2 * inv_r;
  return 3. * r_x * r_z * inv_r5;
}

static inline double D_011(double r_x, double r_y, double r_z, double inv_r) {  // d2/dydz phi(r)

  double inv_r2 = inv_r * inv_r;
  double inv_r5 = inv_r2 * inv_r2 * inv_r;
  return 3. * r_y * r_z * inv_r5;
}



/* 3-rd order */
static inline double D_300(double r_x, double r_y, double r_z, double inv_r) {  // d3/dxdxdx phi(r)

  double inv_r2 = inv_r * inv_r;
  double inv_r5 = inv_r2 * inv_r2 * inv_r;
  double inv_r7 = inv_r5 * inv_r2;
  return -15. * r_x * r_x * r_x * inv_r7 + 9. * r_x * inv_r5;
}

static inline double D_030(double r_x, double r_y, double r_z, double inv_r) {  // d3/dydydy phi(r)

  double inv_r2 = inv_r * inv_r;
  double inv_r5 = inv_r2 * inv_r2 * inv_r;
  double inv_r7 = inv_r5 * inv_r2;
  return -15. * r_y * r_y * r_y * inv_r7 + 9. * r_y * inv_r5;
}

static inline double D_003(double r_x, double r_y, double r_z, double inv_r) {  // d3/dzdzdz phi(r)

  double inv_r2 = inv_r * inv_r;
  double inv_r5 = inv_r2 * inv_r2 * inv_r;
  double inv_r7 = inv_r5 * inv_r2;
  return -15. * r_z * r_z * r_z * inv_r7 + 9. * r_z * inv_r5;
}

static inline double D_210(double r_x, double r_y, double r_z, double inv_r) {  // d3/dxdxdy phi(r)

  double inv_r2 = inv_r * inv_r;
  double inv_r5 = inv_r2 * inv_r2 * inv_r;
  double inv_r7 = inv_r5 * inv_r2;
  return -15. * r_x * r_x * r_y * inv_r7 + 3. * r_y * inv_r5;
}

static inline double D_201(double r_x, double r_y, double r_z, double inv_r) {  // d3/dxdxdz phi(r)

  double inv_r2 = inv_r * inv_r;
  double inv_r5 = inv_r2 * inv_r2 * inv_r;
  double inv_r7 = inv_r5 * inv_r2;
  return -15. * r_x * r_x * r_z * inv_r7 + 3. * r_z * inv_r5;
}

static inline double D_120(double r_x, double r_y, double r_z, double inv_r) {  // d3/dxdydy phi(r)

  double inv_r2 = inv_r * inv_r;
  double inv_r5 = inv_r2 * inv_r2 * inv_r;
  double inv_r7 = inv_r5 * inv_r2;
  return -15. * r_x * r_y * r_y * inv_r7 + 3. * r_x * inv_r5;
}

static inline double D_021(double r_x, double r_y, double r_z, double inv_r) {  // d3/dydydz phi(r)

  double inv_r2 = inv_r * inv_r;
  double inv_r5 = inv_r2 * inv_r2 * inv_r;
  double inv_r7 = inv_r5 * inv_r2;
  return -15. * r_y * r_y * r_z * inv_r7 + 3. * r_z * inv_r5;
}

static inline double D_102(double r_x, double r_y, double r_z, double inv_r) {  // d3/dxdzdz phi(r)

  double inv_r2 = inv_r * inv_r;
  double inv_r5 = inv_r2 * inv_r2 * inv_r;
  double inv_r7 = inv_r5 * inv_r2;
  return -15. * r_x * r_z * r_z * inv_r7 + 3. * r_x * inv_r5;
}

static inline double D_012(double r_x, double r_y, double r_z, double inv_r) {  // d3/dydzdz phi(r)

  double inv_r2 = inv_r * inv_r;
  double inv_r5 = inv_r2 * inv_r2 * inv_r;
  double inv_r7 = inv_r5 * inv_r2;
  return -15. * r_y * r_z * r_z * inv_r7 + 3. * r_y * inv_r5;
}


static inline double D_111(double r_x, double r_y, double r_z, double inv_r) {  // d3/dxdydz phi(r)

  double inv_r3 = inv_r * inv_r * inv_r;
  double inv_r7 = inv_r3 * inv_r3 * inv_r;
  return -15. * r_x * r_z * r_z * inv_r7;
}


/************************************************
 * Field tensor calculation for the accelerations
 ************************************************/


static inline void computeFieldTensors_x(double r_x, double r_y, double r_z, double inv_r, struct multipole A, struct fieldTensor* B)
{

  B->F_000 += A.M_000 * D_100(r_x, r_y, r_z, inv_r);

#ifdef QUADRUPOLES
  B->F_000 += A.M_200 * D_300(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_020 * D_120(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_002 * D_102(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_110 * D_210(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_011 * D_111(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_101 * D_201(r_x, r_y, r_z, inv_r);
#endif
  
  B->F_100 += A.M_000 * D_200(r_x, r_y, r_z, inv_r);
  B->F_010 += A.M_000 * D_110(r_x, r_y, r_z, inv_r);
  B->F_001 += A.M_000 * D_101(r_x, r_y, r_z, inv_r);
  
#ifdef QUADRUPOLES
  B->F_200 += A.M_000 * D_300(r_x, r_y, r_z, inv_r);
  B->F_020 += A.M_000 * D_120(r_x, r_y, r_z, inv_r);
  B->F_002 += A.M_000 * D_102(r_x, r_y, r_z, inv_r);
  B->F_110 += A.M_000 * D_210(r_x, r_y, r_z, inv_r);
  B->F_101 += A.M_000 * D_201(r_x, r_y, r_z, inv_r);
  B->F_011 += A.M_000 * D_111(r_x, r_y, r_z, inv_r);
#endif

}

static inline void computeFieldTensors_y(double r_x, double r_y, double r_z, double inv_r, struct multipole A, struct fieldTensor* B)
{

  B->F_000 += A.M_000 * D_010(r_x, r_y, r_z, inv_r);

#ifdef QUADRUPOLES
  B->F_000 += A.M_200 * D_210(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_020 * D_030(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_002 * D_012(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_110 * D_120(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_011 * D_021(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_101 * D_111(r_x, r_y, r_z, inv_r);
#endif
  
  B->F_100 += A.M_000 * D_110(r_x, r_y, r_z, inv_r);
  B->F_010 += A.M_000 * D_020(r_x, r_y, r_z, inv_r);
  B->F_001 += A.M_000 * D_011(r_x, r_y, r_z, inv_r);
  
#ifdef QUADRUPOLES
  B->F_200 += A.M_000 * D_210(r_x, r_y, r_z, inv_r);
  B->F_020 += A.M_000 * D_030(r_x, r_y, r_z, inv_r);
  B->F_002 += A.M_000 * D_012(r_x, r_y, r_z, inv_r);
  B->F_110 += A.M_000 * D_120(r_x, r_y, r_z, inv_r);
  B->F_101 += A.M_000 * D_111(r_x, r_y, r_z, inv_r);
  B->F_011 += A.M_000 * D_021(r_x, r_y, r_z, inv_r);
#endif

}

static inline void computeFieldTensors_z(double r_x, double r_y, double r_z, double inv_r, struct multipole A, struct fieldTensor* B)
{

  B->F_000 += A.M_000 * D_001(r_x, r_y, r_z, inv_r);

#ifdef QUADRUPOLES
  B->F_000 += A.M_200 * D_201(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_020 * D_021(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_002 * D_003(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_110 * D_111(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_011 * D_012(r_x, r_y, r_z, inv_r);
  B->F_000 += A.M_101 * D_102(r_x, r_y, r_z, inv_r);
#endif
  
  B->F_100 += A.M_000 * D_101(r_x, r_y, r_z, inv_r);
  B->F_010 += A.M_000 * D_011(r_x, r_y, r_z, inv_r);
  B->F_001 += A.M_000 * D_002(r_x, r_y, r_z, inv_r);
  
#ifdef QUADRUPOLES
  B->F_200 += A.M_000 * D_201(r_x, r_y, r_z, inv_r);
  B->F_020 += A.M_000 * D_021(r_x, r_y, r_z, inv_r);
  B->F_002 += A.M_000 * D_003(r_x, r_y, r_z, inv_r);
  B->F_110 += A.M_000 * D_111(r_x, r_y, r_z, inv_r);
  B->F_101 += A.M_000 * D_102(r_x, r_y, r_z, inv_r);
  B->F_011 += A.M_000 * D_012(r_x, r_y, r_z, inv_r);
#endif

}



static inline void shiftAndAddTensor(struct fieldTensor A, struct fieldTensor* B, double dx, double dy, double dz)
{
  B->F_000 += A.F_000;
  B->F_000 += dx * A.F_100;
  B->F_000 += dy * A.F_010; 
  B->F_000 += dz * A.F_001;

#ifdef QUADRUPOLES
  B->F_000 += 0.5f * dx * dx * A.F_200;
  B->F_000 += 0.5f * dy * dy * A.F_020;
  B->F_000 += 0.5f * dz * dz * A.F_002;
  B->F_000 += 0.5f * dx * dy * A.F_110;
  B->F_000 += 0.5f * dx * dz * A.F_101;
  B->F_000 += 0.5f * dy * dz * A.F_011;
#endif

  B->F_100 += A.F_100;
  B->F_010 += A.F_010; 
  B->F_001 += A.F_001; 

#ifdef QUADRUPOLES
  B->F_100 += dx * A.F_200 + dy * A.F_110 + dz * A.F_101;
  B->F_010 += dx * A.F_110 + dy * A.F_020 + dz * A.F_011;
  B->F_001 += dx * A.F_101 + dy * A.F_011 + dz * A.F_002;

  B->F_200 += A.F_200;
  B->F_020 += A.F_020;
  B->F_002 += A.F_002;
  B->F_110 += A.F_110;
  B->F_101 += A.F_101;
  B->F_011 += A.F_011;
#endif
  
}




static inline float applyFieldAcceleration(struct fieldTensor B, double dx, double dy, double dz)
{
  float a = 0.f;

  a += B.F_000;

  //a += dx * B.F_100;
  //a += dy * B.F_010;
  //a += dz * B.F_001;

#ifdef QUADRUPOLES
  a += 0.5f * dx * dx * B.F_200;
  a += 0.5f * dy * dy * B.F_020;
  a += 0.5f * dz * dz * B.F_002;
  a += 0.5f * dx * dy * B.F_110;
  a += 0.5f * dx * dz * B.F_101;
  a += 0.5f * dy * dz * B.F_011;
#endif

  return a;
}