barnes/grav.C at main · ParCoreLab/barnes · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
/*************************************************************************/
/*                                                                       */
/*  Copyright (c) 1994 Stanford University                               */
/*                                                                       */
/*  All rights reserved.                                                 */
/*                                                                       */
/*  Permission is given to use, copy, and modify this software for any   */
/*  non-commercial purpose as long as this copyright notice is not       */
/*  removed.  All other uses, including redistribution in whole or in    */
/*  part, are forbidden without prior written permission.                */
/*                                                                       */
/*  This software is provided with absolutely no warranty and no         */
/*  support.                                                             */
/*                                                                       */
/*************************************************************************/

/*
 * GRAV.C:
 */

EXTERN_ENV
#define global extern

#include "stdinc.h"

#if 0
/*
 * HACKGRAV: evaluate grav field at a given particle.
 */

void hackgrav(bodyptr p, long ProcessId)
{
   Local[ProcessId].pskip = p;
   SETV(Local[ProcessId].pos0, Pos(p));
   Local[ProcessId].phi0 = 0.0;
   CLRV(Local[ProcessId].acc0);
   Local[ProcessId].myn2bterm = 0;
   Local[ProcessId].mynbcterm = 0;
   Local[ProcessId].skipself = FALSE;
   hackwalk(ProcessId);
   Phi(p) = Local[ProcessId].phi0;
   SETV(Acc(p), Local[ProcessId].acc0);
#ifdef QUADPOLE
   Cost(p) = Local[ProcessId].myn2bterm + NDIM * Local[ProcessId].mynbcterm;
#else
   Cost(p) = Local[ProcessId].myn2bterm + Local[ProcessId].mynbcterm;
#endif
}
#endif


/*
 * GRAVSUB: compute a single body-body or body-cell longeraction.
 */

void gravsub(register nodeptr p, long ProcessId)
{
    real drabs, phii, mor3;
    vector ai;

    if (p != Local[ProcessId].pmem) {
        SUBV(Local[ProcessId].dr, Pos(p), Local[ProcessId].pos0);
        DOTVP(Local[ProcessId].drsq, Local[ProcessId].dr, Local[ProcessId].dr);
    }

    Local[ProcessId].drsq += epssq;
    drabs = sqrt((double) Local[ProcessId].drsq);
    phii = Mass(p) / drabs;
    Local[ProcessId].phi0 -= phii;
    mor3 = phii / Local[ProcessId].drsq;
    MULVS(ai, Local[ProcessId].dr, mor3);
    ADDV(Local[ProcessId].acc0, Local[ProcessId].acc0, ai);
    if(Type(p) != BODY) {                  /* a body-cell/leaf interaction? */
       Local[ProcessId].mynbcterm++;
#ifdef QUADPOLE
       dr5inv = 1.0/(Local[ProcessId].drsq * Local[ProcessId].drsq * drabs);
       MULMV(quaddr, Quad(p), Local[ProcessId].dr);
       DOTVP(drquaddr, Local[ProcessId].dr, quaddr);
       phiquad = -0.5 * dr5inv * drquaddr;
       Local[ProcessId].phi0 += phiquad;
       phiquad = 5.0 * phiquad / Local[ProcessId].drsq;
       MULVS(ai, Local[ProcessId].dr, phiquad);
       SUBV(Local[ProcessId].acc0, Local[ProcessId].acc0, ai);
       MULVS(quaddr, quaddr, dr5inv);
       SUBV(Local[ProcessId].acc0, Local[ProcessId].acc0, quaddr);
#endif
    }
    else {                                      /* a body-body interaction  */
       Local[ProcessId].myn2bterm++;
    }
}

#if 0
/*
 * HACKWALK: walk the tree opening cells too close to a given point.
 */

void hackwalk(long ProcessId)
{
    walksub(reinterpret_cast <nodeptr>(Global->G_root), Global->rsize * Global->rsize, ProcessId);
}
#endif

#if 0
/*
 * WALKSUB: recursive routine to do hackwalk operation.
 */

void walksub(nodeptr n, real dsq, long ProcessId)
{
   nodeptr* nn;
   leafptr l;
   bodyptr p;
   long i;

   if (subdivp(n, dsq, ProcessId)) {
      if (Type(n) == CELL) {
	 for (nn = Subp(n); nn < Subp(n) + NSUB; nn++) {
	    if (*nn != NULL) {
	       walksub(*nn, dsq / 4.0, ProcessId);
	    }
	 }
      }
      else {
	 l = (leafptr) n;
	 for (i = 0; i < l->num_bodies; i++) {
	    p = Bodyp(l)[i];
	    if (p != Local[ProcessId].pskip) {
	       gravsub(reinterpret_cast <nodeptr>(p), ProcessId);
	    }
	    else {
	       Local[ProcessId].skipself = TRUE;
	    }
	 }
      }
   }
   else {
      gravsub(n, ProcessId);
   }
}

/*
 * SUBDIVP: decide if a node should be opened.
 * Side effects: sets  pmem,dr, and drsq.
 */

bool subdivp(register nodeptr p, real dsq, long ProcessId)
{
   SUBV(Local[ProcessId].dr, Pos(p), Local[ProcessId].pos0);
   DOTVP(Local[ProcessId].drsq, Local[ProcessId].dr, Local[ProcessId].dr);
   Local[ProcessId].pmem = p;
   return (tolsq * Local[ProcessId].drsq < dsq);
}
#endif