potential.c File Reference

contains functions that compute the potential More...

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include "global.h"
#include "inline.h"

Include dependency graph for potential.c:

Go to the source code of this file.

Functions
double	computeCCPotential (int p1, int p2, int mode)
void	potentialTraverseSubtree (int p, struct bht_node *no)
void	computeRemoteCellsPotential (int rp, struct bht_node *no)
void	computePotentialField ()
void	computePotential ()

Detailed Description

contains functions that compute the potential

Definition in file potential.c.

Function Documentation

double computeCCPotential	(	int	p1,
		int	p2,
		int	mode
	)

This function computes potential for two neighbour cells. mode=0 - computations for local cells mode=1 - computations for remote cells

Definition at line 40 of file potential.c.

References cells, csize, cellData::density, h, h2, h3, statisticsData::mindist, recvData, sdim, cellData::size, partData::size, statistics, x, cellData::x, partData::x, cellData::y, partData::y, cellData::young, partData::young, cellData::z, and partData::z.

{
  double pot;
  double dist;
  double size;
  double g;
  double x, y, z;
  double xc;
  double D;
  double poisson = 0.33;
  double young;

  if (mode == 0 && p1 == p2)
    return 0.0;

  if (mode == 0) {
    x = cells[p1].x;
    y = cells[p1].y;
    z = cells[p1].z;
    size = cells[p1].size;
    young = cells[p1].young;
  }
  if (mode == 1) {
    x = recvData[p1].x;
    y = recvData[p1].y;
    z = recvData[p1].z;
    size = recvData[p1].size;
    young = recvData[p1].young;
  }

  /* compute the distance between two cells */
  if (sdim == 2)
    dist =
    sqrt((x - cells[p2].x) * (x - cells[p2].x) +
         (y - cells[p2].y) * (y - cells[p2].y));
  if (sdim == 3)
    dist =
    sqrt((x - cells[p2].x) * (x - cells[p2].x) +
         (y - cells[p2].y) * (y - cells[p2].y) + (z -
                              cells[p2].z) * (z -
                                      cells
                                      [p2].
                                      z));

  if (statistics.mindist > dist)
    statistics.mindist = dist;

  if (dist <= h) {

    double r01, r02;
    double area;
    double sc;

    if (sdim == 2)
      sc = h2;
    if (sdim == 3)
      sc = h3;

    if (mode == 0) {
      cells[p1].density +=
      sc * (cells[p2].size / csize) * sph_kernel(dist);
    }
    if (mode == 1) {
      cells[p2].density += sc * (size / csize) * sph_kernel(dist);
    }

    xc = size + cells[p2].size - dist;

    if (xc <= 0.0)
      return 0.0;

    D = 0.75 * ((1.0 - poisson * poisson) / young +
        (1.0 - poisson * poisson / cells[p2].young));

    /* adhesion */
    r01 =
    (size * size - cells[p2].size * cells[p2].size +
     dist * dist) / (2 * dist);
    r02 = dist - r01;

    area =
    M_PI *
    ((size * size * (size - r01) -
      (size * size * size - r01 * r01 * r01) / 3) +
     (cells[p2].size * cells[p2].size * (cells[p2].size - r02) -
      (cells[p2].size * cells[p2].size * cells[p2].size -
       r02 * r02 * r02) / 3));

    /* compute potential */
    pot =
    (2.0 * pow(xc, 5 / 2) / (5.0 * D)) * sqrt((size * cells[p2].size) /
                          (size +
                           cells[p2].size)) +
    area * 0.1;

    return pot;

  } else
    return 0.0;

}

Here is the caller graph for this function:

void computePotential

(

)

This is a driving routine for the potential computations.

Definition at line 287 of file potential.c.

References computePotentialField(), MPIrank, statOutStep, and step.

{
#ifdef DEBUG
  if (MPIrank == 0 && !(step % statOutStep)) {
    printf(" Potential computations...");
    fflush(stdout);
  }
#endif

  computePotentialField();

#ifdef DEBUG
  if (MPIrank == 0 && !(step % statOutStep)) {
    printf("done\n");
    fflush(stdout);
  }
#endif
}

Here is the call graph for this function:

Here is the caller graph for this function:

void computePotentialField

(

)

This function computes potential for all local cells.

Definition at line 209 of file potential.c.

References bht_root, cells, cellsExchangeInit(), cellsExchangeWait(), bht_node::child, computeRemoteCellsPotential(), cellData::density, bht_node::father, h, leafs, lnc, statisticsData::maxsize, statisticsData::maxvel, statisticsData::mindist, statisticsData::minsize, statisticsData::minvel, nc, numImp, nx, bht_node::partnum, potentialTraverseSubtree(), s, statistics, tnc, cellData::v, and x.

{
  int p, rp;
  int s;
  struct bht_node *no;

  if (nc > 1)
    statistics.mindist = nx;    /* minimal distance is set to the box size */
  else
    statistics.mindist = 0.0;
  statistics.minvel = DBL_MAX;  /* minimal velocity is set to DBL_MAX */
  statistics.maxvel = 0;    /* maximal velocity is set to zero */
  statistics.minsize = DBL_MAX; /* minimal size is set to DBL_MAX */
  statistics.maxsize = 0;   /* maximal size is set to zero */

  if (lnc == 0)
    return;

  /* initiate asynchronous data transfers between processors */
  cellsExchangeInit();

  /* loop over all local cells */
  //#pragma omp parallel for private(p,s,no) schedule(dynamic)
#pragma omp parallel for private(p,s,no)
  for (p = 0; p < lnc; p++) {

    cells[p].density = 0.0;
    cells[p].v = 0.0;

    /* set the pointer to the node that owns the cell */
    no = leafs[p];

    /* traverse the tree (bottom-up) */
    while (1) {

      if (no->father == NULL)
    break;

      /* loop over siblings */
      for (s = 0; s < tnc; s++) {
    /* "Son of my father but not my brother. It's me." Hairdresser */
    if ((no->father->child[s] == NULL) || (no->father->child[s] == no)
        || (no->father->child[s]->partnum == -1))
      continue;
    if (cellIntersectNode
        (cells[p].x, cells[p].y, cells[p].z, cells[p].h,
         no->father->child[s])) {
      /* traverse subtree of the sibling (up-bottom) */
      potentialTraverseSubtree(p, no->father->child[s]);
    }
      }

      no = no->father;
      if (no->father == NULL)
    break;          /* this is the root node */
      if (cellInsideNode(p, no))
    break;

    }

  }

  /* wait for data transfers to finish */
  cellsExchangeWait();

  /* loop over remote cells */
  for (rp = 0; rp < numImp; rp++) {
    /* start with the root node */
    no = bht_root;
    /* traverse the tree and compute potential (starting from the root) */
    computeRemoteCellsPotential(rp, no);

  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void computeRemoteCellsPotential	(	int	rp,
		struct bht_node *	no
	)

This function traverse the tree (up-bottom) and computes potential for remote cells.

Definition at line 176 of file potential.c.

References cells, bht_node::child, computeCCPotential(), h, bht_node::leaf, lnc, bht_node::partnum, recvData, s, tnc, cellData::v, and x.

{
  int s;

  if (lnc == 0)
    return;
  if (no->partnum == -1)
    return;

  /* compute potential if it is a leaf node */
  if (no->leaf == 1) {
    if (lnc > 0)
      cells[no->partnum].v += computeCCPotential(rp, no->partnum, 1);
    return;
  }

  /* loop over all siblings */
  for (s = 0; s < tnc; s++) {
    /* continue if there is no child */
    if (no->child[s] == NULL || no->child[s]->partnum == -1)
      continue;
    /* proceed if the child node and particle intersect */
    if (!cellIntersectNode
    (recvData[rp].x, recvData[rp].y, recvData[rp].z, recvData[rp].h,
     no->child[s]))
      continue;
    computeRemoteCellsPotential(rp, no->child[s]);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void potentialTraverseSubtree	(	int	p,
		struct bht_node *	no
	)

This function traverses the subtree under given node looking for neighbour cells.

Definition at line 146 of file potential.c.

References cells, bht_node::child, computeCCPotential(), h, bht_node::leaf, lnc, bht_node::partnum, s, tnc, cellData::v, and x.

{
  int s;

  if (no->partnum == -1)
    return;

  /* if it is a leaf node we do not proceed to the siblings */
  if (no->leaf == 1) {
    if (lnc > 0)
      cells[p].v += computeCCPotential(p, no->partnum, 0);
    return;
  }

  /* loop over siblings */
  for (s = 0; s < tnc; s++) {
    /* if it is an empty child proceed to the next one */
    if ((no->child[s] == NULL) || (no->child[s]->partnum == -1))
      continue;
    /* if there is no intersection with this child proceed to next one */
    if (!cellIntersectNode
    (cells[p].x, cells[p].y, cells[p].z, cells[p].h, no->child[s]))
      continue;
    potentialTraverseSubtree(p, no->child[s]);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function: