Timothy  0.9
Tissue Modelling Framework
 All Data Structures Files Functions Variables Typedefs Macros
gradient.c
Go to the documentation of this file.
1 /* **************************************************************************
2  * This file is part of Timothy
3  *
4  * Copyright (c) 2014/15 Maciej Cytowski
5  * Copyright (c) 2014/15 ICM, University of Warsaw, Poland
6  *
7  * This program is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program; if not, write to the Free Software
19  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20  *
21  * *************************************************************************/
22 
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <math.h>
26 #include <float.h>
27 
28 #include "global.h"
29 #include "inline.h"
30 
40 void computeCCgradient(int p1, int p2, int mode)
41 {
42 
43  double pot;
44  double dist;
45  double g;
46  double grad[3];
47  double m = 1.0; /* we assume that cells' mass is always 1.0 */
48  double sc;
49  double x1, x2, y1, y2, z1, z2;
50  double v, density, size;
51  double r;
52 
53  if (p1 == p2 && mode == 0)
54  return;
55 
56  if (mode == 0) {
57  x1 = cells[p1].x;
58  x2 = cells[p2].x;
59  y1 = cells[p1].y;
60  y2 = cells[p2].y;
61  z1 = cells[p1].z;
62  z2 = cells[p2].z;
63  v = cells[p2].v;
64  density = cells[p2].density;
65  size = cells[p2].size;
66  }
67  if (mode == 1) {
68  x1 = recvData[p1].x;
69  x2 = cells[p2].x;
70  y1 = recvData[p1].y;
71  y2 = cells[p2].y;
72  z1 = recvData[p1].z;
73  z2 = cells[p2].z;
74  v = recvDensPotData[p1].v;
75  density = recvDensPotData[p1].density;
76  size = recvData[p1].size;
77  }
78 
79  if (sdim == 2)
80  r = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
81  if (sdim == 3)
82  r = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) +
83  (z1 - z2) * (z1 - z2));
84 
85  grad[0] = 0.0;
86  grad[1] = 0.0;
87  grad[2] = 0.0;
88 
89  /* compute the gradient of SPH kernel function */
90  sph_kernel_gradient(p1, p2, grad, mode, r);
91 
92  if (density == 0.0)
93  sc = 0.0;
94  else {
95  m = size / csize;
96  sc = (m / density) * v;
97  }
98 
99  /* update forces */
100  if (mode == 0) {
101  velocity[p1].x += sc * grad[0];
102  velocity[p1].y += sc * grad[1];
103  velocity[p1].z += sc * grad[2];
104  } else {
105  velocity[p2].x -= sc * grad[0];
106  velocity[p2].y -= sc * grad[1];
107  velocity[p2].z -= sc * grad[2];
108  }
109 
110 }
111 
115 void gradientTraverseSubtree(int p, struct bht_node *no)
116 {
117 
118  int s;
119 
120  if (no->partnum == -1)
121  return;
122 
123  /* if it is a leaf node we do not proceed to the siblings */
124  if (no->leaf == 1) {
125  if (lnc > 0)
126  computeCCgradient(p, no->partnum, 0);
127  return;
128  }
129 
130  /* loop over siblings */
131  for (s = 0; s < tnc; s++) {
132  /* if it is an empty child proceed to the next one */
133  if (no->child[s] == NULL || no->child[s]->partnum == -1)
134  continue;
135  /* if there is no intersection with this child proceed to next one */
136  if (!cellIntersectNode
137  (cells[p].x, cells[p].y, cells[p].z, cells[p].h, no->child[s]))
138  continue;
139  /* if it is an leaf child compute the potential */
140  if (no->child[s]->leaf == 1) {
141  computeCCgradient(p, no->child[s]->partnum, 0);
142  } else {
143  gradientTraverseSubtree(p, no->child[s]);
144  }
145  }
146 }
147 
152 {
153 
154  int s;
155 
156  if (lnc == 0)
157  return;
158  if (no->partnum == -1)
159  return;
160 
161  /* compute potential if it is a leaf node */
162  if (no->leaf == 1) {
163  computeCCgradient(rp, no->partnum, 1);
164  return;
165  }
166 
167  /* loop over all siblings */
168  for (s = 0; s < tnc; s++) {
169  /* continue if there is no child */
170  if (no->child[s] == NULL || no->child[s]->partnum == -1)
171  continue;
172  /* proceed if the child node and particle intersect */
173  if (cellIntersectNode
174  (recvData[rp].x, recvData[rp].y, recvData[rp].z, recvData[rp].h,
175  no->child[s]))
177  }
178 }
179 
184 {
185  int p, s;
186  struct bht_node *no;
187  double dvel;
188  double sf;
189 #ifdef DEBUG
190  if (MPIrank == 0 && !(step % statOutStep)) {
191  printf(" Potential gradient...");
192  fflush(stdout);
193  }
194 #endif
196 
197  /* loop over all local cells */
198 #pragma omp parallel for private(p,s,no) //schedule(dynamic)
199  for (p = 0; p < lnc; p++) {
200  velocity[p].x = 0.0;
201  velocity[p].y = 0.0;
202  velocity[p].z = 0.0;
203  /* set the pointer to the node that owns the cell */
204  no = leafs[p];
205  /* traverse the tree (bottom-up) */
206  while (1) {
207  if (no->father == NULL)
208  break;
209  /* loop over siblings */
210  for (s = 0; s < tnc; s++) {
211  /* "Son of my father but not my brother. It's me." Hairdresser */
212  if ((no->father->child[s] == NULL) || (no->father->child[s] == no)
213  || (no->father->child[s]->partnum == -1))
214  continue;
215  if (cellIntersectNode
216  (cells[p].x, cells[p].y, cells[p].z, cells[p].h,
217  no->father->child[s])) {
218  /* traverse subtree of the sibling (up-bottom) */
220  }
221 
222  }
223  no = no->father;
224  if (no->father == NULL)
225  break; /* this is the root node */
226  if (cellInsideNode(p, no))
227  break;
228  }
229  }
230 
232 
233  /* loop over imported cells */
234  for (p = 0; p < numImp; p++) {
235  /* start with the root node */
236  no = bht_root;
237  /* traverse the tree starting from the root and compute potential */
239 
240  }
241 
242  for (p = 0; p < lnc; p++) {
243  dvel =
244  sqrt(velocity[p].x * velocity[p].x +
245  velocity[p].y * velocity[p].y +
246  velocity[p].z * velocity[p].z);
247  if (dvel < statistics.minvel)
248  statistics.minvel = dvel;
249  if (dvel > statistics.maxvel)
250  statistics.maxvel = dvel;
251  if (cells[p].size < statistics.minsize)
253  if (cells[p].size > statistics.maxsize)
255  }
256  /* this should be removed soon (do we really need to reduceall here?) */
257  MPI_Allreduce(&statistics.minvel, &globalMinVel, 1, MPI_DOUBLE, MPI_MIN,
258  MPI_COMM_WORLD);
259  MPI_Allreduce(&statistics.maxvel, &globalMaxVel, 1, MPI_DOUBLE, MPI_MAX,
260  MPI_COMM_WORLD);
261 
262  if (globalMaxVel == 0.0)
263  sf = 0.0;
264  else
266 
267  statistics.minvel = DBL_MAX; /* minimal velocity is set to DBL_MAX */
268  statistics.maxvel = 0; /* maximal velocity is set to zero */
269 
270  for (p = 0; p < lnc; p++) {
271  velocity[p].x *= sf;
272  velocity[p].y *= sf;
273  velocity[p].z *= sf;
274  dvel =
275  sqrt(velocity[p].x * velocity[p].x +
276  velocity[p].y * velocity[p].y +
277  velocity[p].z * velocity[p].z);
278  if (dvel < statistics.minvel)
279  statistics.minvel = dvel;
280  if (dvel > statistics.maxvel)
281  statistics.maxvel = dvel;
282  }
283 #ifdef DEBUG
284  if (MPIrank == 0 && !(step % statOutStep)) {
285  printf("done\n");
286  fflush(stdout);
287  }
288 #endif
289 }
double density
Definition: global.h:75
int statOutStep
Definition: global.h:214
double z
Definition: global.h:225
float secondsPerStep
Definition: global.h:247
void gradientTraverseSubtree(int p, struct bht_node *no)
Definition: gradient.c:115
int MPIrank
Definition: global.h:134
#define lnc
Definition: global.h:102
float s
Definition: global.h:186
HYPRE_SStructVector x
Definition: tempf.c:41
unsigned char leaf
Definition: global.h:284
double z
Definition: global.h:119
contains various Timothy inline functions
void densPotExchangeInit()
Definition: comm.c:257
double maxsize
Definition: global.h:258
double size
Definition: global.h:72
struct partData * recvData
Definition: global.h:149
contains the most important global variables, arrays and defines
struct doubleVector3d * velocity
Definition: global.h:84
float m
Definition: global.h:188
void computeCCgradient(int p1, int p2, int mode)
Definition: gradient.c:40
void densPotExchangeWait()
Definition: comm.c:301
int sdim
Definition: global.h:160
double minsize
Definition: global.h:257
int tnc
Definition: global.h:295
int step
Definition: global.h:173
double y
Definition: global.h:224
struct bht_node * child[8]
Definition: global.h:283
struct statisticsData statistics
Definition: global.h:268
double x
Definition: global.h:69
double x
Definition: global.h:223
int partnum
Definition: global.h:281
double csize
Definition: global.h:197
struct bht_node * bht_root
Definition: global.h:289
struct bht_node ** leafs
Definition: global.h:290
int numImp
Definition: global.h:154
struct cellData * cells
Definition: global.h:82
struct bht_node * father
Definition: global.h:282
void computeRemoteCellsPotentialGradient(int rp, struct bht_node *no)
Definition: gradient.c:151
double v
Definition: global.h:127
float v
Definition: global.h:189
double size
Definition: global.h:121
double maxvel
Definition: global.h:260
double z
Definition: global.h:71
struct densPotData * recvDensPotData
Definition: global.h:151
double y
Definition: global.h:70
double globalMinVel
Definition: global.h:270
double globalMaxVel
Definition: global.h:271
double x
Definition: global.h:117
double v
Definition: global.h:74
double y
Definition: global.h:118
float maxSpeedInUnits
Definition: global.h:177
double density
Definition: global.h:128
double minvel
Definition: global.h:261
double h
Definition: global.h:200
int computePotentialGradient()
Definition: gradient.c:183