Skip to content
Snippets Groups Projects
Commit f5804e06 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Added face calculations and unit test for 3D Voronoi algorithm. Started...

Added face calculations and unit test for 3D Voronoi algorithm. Started testing 3D moving mesh, does not work.
parent 09e08a22
No related branches found
No related tags found
1 merge request!3211D and 2D moving mesh algorithm
......@@ -37,9 +37,9 @@
#define const_max_u_change 0.1f
/* Dimensionality of the problem */
//#define HYDRO_DIMENSION_3D
#define HYDRO_DIMENSION_3D
//#define HYDRO_DIMENSION_2D
#define HYDRO_DIMENSION_1D
//#define HYDRO_DIMENSION_1D
/* Hydrodynamical adiabatic index. */
#define HYDRO_GAMMA_5_3
......
......@@ -86,15 +86,14 @@ __attribute__((always_inline)) INLINE void hydro_gradients_single_quantity(
__attribute__((always_inline)) INLINE static void hydro_gradients_collect(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
float A[3], midpoint[3];
float A, midpoint[3];
if (!voronoi_get_face(&pi->cell, pj->id, A, midpoint)) {
A = voronoi_get_face(&pi->cell, pj->id, midpoint);
if (!A) {
/* particle is not a cell neighbour: do nothing */
return;
}
float Anorm = sqrtf(A[0] * A[0] + A[1] * A[1] + A[2] * A[2]);
float c[3];
/* midpoint is relative w.r.t. pi->x, as is dx */
/* c is supposed to be the vector pointing from the midpoint of pi and pj to
......@@ -108,15 +107,15 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
float r = sqrtf(r2);
hydro_gradients_single_quantity(pi->primitives.rho, pj->primitives.rho, c, dx,
r, Anorm, pi->primitives.gradients.rho);
r, A, pi->primitives.gradients.rho);
hydro_gradients_single_quantity(pi->primitives.v[0], pj->primitives.v[0], c,
dx, r, Anorm, pi->primitives.gradients.v[0]);
dx, r, A, pi->primitives.gradients.v[0]);
hydro_gradients_single_quantity(pi->primitives.v[1], pj->primitives.v[1], c,
dx, r, Anorm, pi->primitives.gradients.v[1]);
dx, r, A, pi->primitives.gradients.v[1]);
hydro_gradients_single_quantity(pi->primitives.v[2], pj->primitives.v[2], c,
dx, r, Anorm, pi->primitives.gradients.v[2]);
dx, r, A, pi->primitives.gradients.v[2]);
hydro_gradients_single_quantity(pi->primitives.P, pj->primitives.P, c, dx, r,
Anorm, pi->primitives.gradients.P);
A, pi->primitives.gradients.P);
hydro_slope_limit_cell_collect(pi, pj, r);
......@@ -125,19 +124,15 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_collect(
mindx[1] = -dx[1];
mindx[2] = -dx[2];
hydro_gradients_single_quantity(pj->primitives.rho, pi->primitives.rho, c,
mindx, r, Anorm,
pj->primitives.gradients.rho);
mindx, r, A, pj->primitives.gradients.rho);
hydro_gradients_single_quantity(pj->primitives.v[0], pi->primitives.v[0], c,
mindx, r, Anorm,
pj->primitives.gradients.v[0]);
mindx, r, A, pj->primitives.gradients.v[0]);
hydro_gradients_single_quantity(pj->primitives.v[1], pi->primitives.v[1], c,
mindx, r, Anorm,
pj->primitives.gradients.v[1]);
mindx, r, A, pj->primitives.gradients.v[1]);
hydro_gradients_single_quantity(pj->primitives.v[2], pi->primitives.v[2], c,
mindx, r, Anorm,
pj->primitives.gradients.v[2]);
mindx, r, A, pj->primitives.gradients.v[2]);
hydro_gradients_single_quantity(pj->primitives.P, pi->primitives.P, c, mindx,
r, Anorm, pj->primitives.gradients.P);
r, A, pj->primitives.gradients.P);
hydro_slope_limit_cell_collect(pj, pi, r);
}
......@@ -156,15 +151,14 @@ __attribute__((always_inline)) INLINE static void
hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
struct part *pi, struct part *pj) {
float A[3], midpoint[3];
float A, midpoint[3];
if (!voronoi_get_face(&pi->cell, pj->id, A, midpoint)) {
A = voronoi_get_face(&pi->cell, pj->id, midpoint);
if (!A) {
/* particle is not a cell neighbour: do nothing */
return;
}
float Anorm = sqrtf(A[0] * A[0] + A[1] * A[1] + A[2] * A[2]);
float c[3];
/* midpoint is relative w.r.t. pi->x, as is dx */
/* c is supposed to be the vector pointing from the midpoint of pi and pj to
......@@ -178,15 +172,15 @@ hydro_gradients_nonsym_collect(float r2, float *dx, float hi, float hj,
float r = sqrtf(r2);
hydro_gradients_single_quantity(pi->primitives.rho, pj->primitives.rho, c, dx,
r, Anorm, pi->primitives.gradients.rho);
r, A, pi->primitives.gradients.rho);
hydro_gradients_single_quantity(pi->primitives.v[0], pj->primitives.v[0], c,
dx, r, Anorm, pi->primitives.gradients.v[0]);
dx, r, A, pi->primitives.gradients.v[0]);
hydro_gradients_single_quantity(pi->primitives.v[1], pj->primitives.v[1], c,
dx, r, Anorm, pi->primitives.gradients.v[1]);
dx, r, A, pi->primitives.gradients.v[1]);
hydro_gradients_single_quantity(pi->primitives.v[2], pj->primitives.v[2], c,
dx, r, Anorm, pi->primitives.gradients.v[2]);
dx, r, A, pi->primitives.gradients.v[2]);
hydro_gradients_single_quantity(pi->primitives.P, pj->primitives.P, c, dx, r,
Anorm, pi->primitives.gradients.P);
A, pi->primitives.gradients.P);
hydro_slope_limit_cell_collect(pi, pj, r);
}
......
......@@ -134,8 +134,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float r = sqrtf(r2);
int k;
float A[3];
float Anorm;
float A;
float xij_i[3];
float vmax, dvdotdx;
float vi[3], vj[3], vij[3];
......@@ -143,7 +142,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float dti, dtj, mindt;
float n_unit[3];
if (!voronoi_get_face(&pi->cell, pj->id, A, xij_i)) {
A = voronoi_get_face(&pi->cell, pj->id, xij_i);
if (!A) {
/* this neighbour does not share a face with the cell, return */
return;
}
......@@ -190,11 +190,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
dti = mindt;
dtj = mindt;
Anorm = A[0] * A[0] + A[1] * A[1] + A[2] * A[2];
/* compute the normal vector of the interface */
Anorm = sqrtf(Anorm);
for (k = 0; k < 3; k++) n_unit[k] = A[k] / Anorm;
for (k = 0; k < 3; ++k) {
n_unit[k] = -dx[k] / r;
}
/* Compute interface velocity */
float fac = (vi[0] - vj[0]) * (xij_i[0] + 0.5f * dx[0]) +
......@@ -261,19 +260,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
/* Update conserved variables */
/* eqn. (16) */
pi->conserved.flux.mass -= dti * Anorm * totflux[0];
pi->conserved.flux.momentum[0] -= dti * Anorm * totflux[1];
pi->conserved.flux.momentum[1] -= dti * Anorm * totflux[2];
pi->conserved.flux.momentum[2] -= dti * Anorm * totflux[3];
pi->conserved.flux.energy -= dti * Anorm * totflux[4];
pi->conserved.flux.mass -= dti * A * totflux[0];
pi->conserved.flux.momentum[0] -= dti * A * totflux[1];
pi->conserved.flux.momentum[1] -= dti * A * totflux[2];
pi->conserved.flux.momentum[2] -= dti * A * totflux[3];
pi->conserved.flux.energy -= dti * A * totflux[4];
float ekin = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
pi->primitives.v[1] * pi->primitives.v[1] +
pi->primitives.v[2] * pi->primitives.v[2]);
pi->conserved.flux.energy += dti * Anorm * totflux[1] * pi->primitives.v[0];
pi->conserved.flux.energy += dti * Anorm * totflux[2] * pi->primitives.v[1];
pi->conserved.flux.energy += dti * Anorm * totflux[3] * pi->primitives.v[2];
pi->conserved.flux.energy -= dti * Anorm * totflux[0] * ekin;
pi->conserved.flux.energy += dti * A * totflux[1] * pi->primitives.v[0];
pi->conserved.flux.energy += dti * A * totflux[2] * pi->primitives.v[1];
pi->conserved.flux.energy += dti * A * totflux[3] * pi->primitives.v[2];
pi->conserved.flux.energy -= dti * A * totflux[0] * ekin;
/* here is how it works:
Mode will only be 1 if both particles are ACTIVE and they are in the same
......@@ -288,19 +287,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
==> we update particle j if (MODE IS 1) OR (j IS INACTIVE)
*/
if (mode == 1 || pj->ti_end > pi->ti_end) {
pj->conserved.flux.mass += dtj * Anorm * totflux[0];
pj->conserved.flux.momentum[0] += dtj * Anorm * totflux[1];
pj->conserved.flux.momentum[1] += dtj * Anorm * totflux[2];
pj->conserved.flux.momentum[2] += dtj * Anorm * totflux[3];
pj->conserved.flux.energy += dtj * Anorm * totflux[4];
pj->conserved.flux.mass += dtj * A * totflux[0];
pj->conserved.flux.momentum[0] += dtj * A * totflux[1];
pj->conserved.flux.momentum[1] += dtj * A * totflux[2];
pj->conserved.flux.momentum[2] += dtj * A * totflux[3];
pj->conserved.flux.energy += dtj * A * totflux[4];
ekin = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
pj->primitives.v[1] * pj->primitives.v[1] +
pj->primitives.v[2] * pj->primitives.v[2]);
pj->conserved.flux.energy -= dtj * Anorm * totflux[1] * pj->primitives.v[0];
pj->conserved.flux.energy -= dtj * Anorm * totflux[2] * pj->primitives.v[1];
pj->conserved.flux.energy -= dtj * Anorm * totflux[3] * pj->primitives.v[2];
pj->conserved.flux.energy += dtj * Anorm * totflux[0] * ekin;
pj->conserved.flux.energy -= dtj * A * totflux[1] * pj->primitives.v[0];
pj->conserved.flux.energy -= dtj * A * totflux[2] * pj->primitives.v[1];
pj->conserved.flux.energy -= dtj * A * totflux[3] * pj->primitives.v[2];
pj->conserved.flux.energy += dtj * A * totflux[0] * ekin;
}
}
......
......@@ -128,43 +128,37 @@ __attribute__((always_inline)) INLINE float voronoi_cell_finalize(
* deal with it appropriately.
*
* For this specific case, we simply check if the neighbour is the left or
* right neighbour and set the surface area to the unit vector pointing either
* to the right or to the left. The midpoint is set to the relative position
* vector of the appropriate face.
* right neighbour and set the surface area to 1. The midpoint is set to the
* relative position vector of the appropriate face.
*
* @param cell 1D Voronoi cell.
* @param ngb ID of a particle that is possibly a neighbour of this cell.
* @param A Array to store the oriented surface area in.
* @param midpoint Array to store the relative position of the face in.
* @return 0 if the given neighbour is not a neighbour, 1 otherwise.
* @return 0 if the given neighbour is not a neighbour, surface area 1.0f
* otherwise.
*/
__attribute__((always_inline)) INLINE int voronoi_get_face(
struct voronoi_cell *cell, unsigned long long ngb, float *A,
float *midpoint) {
__attribute__((always_inline)) INLINE float voronoi_get_face(
struct voronoi_cell *cell, unsigned long long ngb, float *midpoint) {
if (ngb != cell->idL && ngb != cell->idR) {
/* this is perfectly possible: we interact with all particles within the
smoothing length, and they do not need to be all neighbours.
If this happens, we return 0, so that the flux method can return */
return 0;
return 0.0f;
}
if (ngb == cell->idL) {
/* Left face */
A[0] = -1.0f;
midpoint[0] = cell->xL;
} else {
/* Right face */
A[0] = 1.0f;
midpoint[0] = cell->xR;
}
/* The other components of A and midpoint are just zero */
A[1] = 0.0f;
A[2] = 0.0f;
/* The other components of midpoint are just zero */
midpoint[1] = 0.0f;
midpoint[2] = 0.0f;
return 1;
return 1.0f;
}
/**
......
......@@ -21,6 +21,9 @@
#define SWIFT_VORONOI3D_ALGORITHM_H
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "error.h"
#include "inline.h"
#include "voronoi3d_cell.h"
......@@ -1219,6 +1222,13 @@ __attribute__((always_inline)) INLINE void voronoi_calculate_cell(
cell->centroid[0] += cell->x[0];
cell->centroid[1] += cell->x[1];
cell->centroid[2] += cell->x[2];
/* Reset the edges: we still need them for the face calculation */
for (i = 0; i < VORONOI3D_MAXNUMEDGE; ++i) {
if (cell->edges[i] < 0) {
cell->edges[i] = -1 - cell->edges[i];
}
}
}
/**
......@@ -1389,24 +1399,36 @@ __attribute__((always_inline)) INLINE float voronoi_cell_finalize(
}
max_radius = sqrtf(max_radius);
return max_radius;
return 2.0f * max_radius;
}
/**
* @brief Get the oriented surface area and midpoint of the face between a
* 3D Voronoi cell and the given neighbour
* @brief Get the surface area and midpoint of the face between a 3D Voronoi
* cell and the given neighbour
*
* @param cell 3D Voronoi cell.
* @param ngb ID of a particle that is possibly a neighbour of this cell->
* @param A Array to store the oriented surface area in.
* @param midpoint Array to store the relative position of the face in.
* @return 0 if the given neighbour is not a neighbour, 1 otherwise.
* @return 0 if the given neighbour is not a neighbour, the surface area of
* the face otherwise.
*/
__attribute__((always_inline)) INLINE int voronoi_get_face(
struct voronoi_cell *cell, unsigned long long ngb, float *A,
float *midpoint) {
__attribute__((always_inline)) INLINE float voronoi_get_face(
struct voronoi_cell *cell, unsigned long long ngb, float *midpoint) {
return 0;
int i = 0;
while (i < cell->nface && cell->ngbs[i] != ngb) {
i++;
}
if (i == cell->nface) {
/* Ngb not found */
return 0.0f;
}
midpoint[0] = cell->face_midpoints[i][0];
midpoint[1] = cell->face_midpoints[i][1];
midpoint[2] = cell->face_midpoints[i][2];
return cell->face_areas[i];
}
/**
......
......@@ -55,18 +55,18 @@ int main() {
}
/* Check face method */
float A[3];
float A;
float midpoint[3];
voronoi_get_face(&cell, 1, A, midpoint);
if (A[0] != -1.0f) {
error("Wrong surface normal returned for left neighbour!");
A = voronoi_get_face(&cell, 1, midpoint);
if (A != 1.0f) {
error("Wrong surface area returned for left neighbour!");
}
if (midpoint[0] != -0.25f) {
error("Wrong midpoint returned for left neighbour!");
}
voronoi_get_face(&cell, 2, A, midpoint);
if (A[0] != 1.0f) {
error("Wrong surface normal returned for right neighbour!");
A = voronoi_get_face(&cell, 2, midpoint);
if (A != 1.0f) {
error("Wrong surface area returned for right neighbour!");
}
if (midpoint[0] != 0.25f) {
error("Wrong midpoint returned for right neighbour!");
......
......@@ -64,12 +64,54 @@ void test_calculate_cell() {
voronoi_initialize(&cell);
/* Calculate the volume and centroid of the large cube. */
voronoi_calculate_cell(&cell);
/* Calculate the faces. */
voronoi_calculate_faces(&cell);
/* Update these values if you ever change to another large cube! */
assert(cell.volume == 27.0f);
assert(cell.centroid[0] = 0.5f);
assert(cell.centroid[1] = 0.5f);
assert(cell.centroid[2] = 0.5f);
/* Check cell neighbours. */
assert(cell.nface == 6);
assert(cell.ngbs[0] == VORONOI3D_BOX_FRONT);
assert(cell.ngbs[1] == VORONOI3D_BOX_LEFT);
assert(cell.ngbs[2] == VORONOI3D_BOX_BOTTOM);
assert(cell.ngbs[3] == VORONOI3D_BOX_TOP);
assert(cell.ngbs[4] == VORONOI3D_BOX_BACK);
assert(cell.ngbs[5] == VORONOI3D_BOX_RIGHT);
/* Check cell faces */
assert(cell.face_areas[0] == 9.0f);
assert(cell.face_midpoints[0][0] == 0.5f);
assert(cell.face_midpoints[0][1] == -1.0f);
assert(cell.face_midpoints[0][2] == 0.5f);
assert(cell.face_areas[1] == 9.0f);
assert(cell.face_midpoints[1][0] == -1.0f);
assert(cell.face_midpoints[1][1] == 0.5f);
assert(cell.face_midpoints[1][2] == 0.5f);
assert(cell.face_areas[2] == 9.0f);
assert(cell.face_midpoints[2][0] == 0.5f);
assert(cell.face_midpoints[2][1] == 0.5f);
assert(cell.face_midpoints[2][2] == -1.0f);
assert(cell.face_areas[3] == 9.0f);
assert(cell.face_midpoints[3][0] == 0.5f);
assert(cell.face_midpoints[3][1] == 0.5f);
assert(cell.face_midpoints[3][2] == 2.0f);
assert(cell.face_areas[4] == 9.0f);
assert(cell.face_midpoints[4][0] == 0.5f);
assert(cell.face_midpoints[4][1] == 2.0f);
assert(cell.face_midpoints[4][2] == 0.5f);
assert(cell.face_areas[5] == 9.0f);
assert(cell.face_midpoints[5][0] == 2.0f);
assert(cell.face_midpoints[5][1] == 0.5f);
assert(cell.face_midpoints[5][2] == 0.5f);
}
int main() {
......@@ -116,14 +158,22 @@ int main() {
error("Wrong centroid: %g %g %g!", cell.centroid[0], cell.centroid[1],
cell.centroid[2]);
}
/* Check neighbour order */
// TODO
/* Check face method */
float A[3];
float midpoint[3];
voronoi_get_face(&cell, 1, A, midpoint);
// TODO
/* Check faces. */
float A, midpoint[3];
A = voronoi_get_face(&cell, 1, midpoint);
if (A) {
if (A != 0.25f) {
error("Wrong surface area: %g!", A);
}
if (fabs(midpoint[0] - 0.25f) > 1.e-5 || fabs(midpoint[1] - 0.5f) > 1.e-5 ||
fabs(midpoint[2] - 0.5f) > 1.e-5) {
error("Wrong face midpoint: %g %g %g!", midpoint[0], midpoint[1],
midpoint[2]);
}
} else {
error("Neighbour not found!");
}
return 0;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment