Add version with stubbed force calculation
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
This commit is contained in:
Rafael Ravedutti
70
src/core/allocate.c
Normal file
70
src/core/allocate.c
Normal file
@@ -0,0 +1,70 @@
|
||||
/*
|
||||
* =======================================================================================
|
||||
*
|
||||
* Author: Jan Eitzinger (je), jan.eitzinger@fau.de
|
||||
* Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
|
||||
*
|
||||
* This file is part of MD-Bench.
|
||||
*
|
||||
* MD-Bench is free software: you can redistribute it and/or modify it
|
||||
* under the terms of the GNU Lesser General Public License as published
|
||||
* by the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
|
||||
* PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
|
||||
* details.
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License along
|
||||
* with MD-Bench. If not, see <https://www.gnu.org/licenses/>.
|
||||
* =======================================================================================
|
||||
*/
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
#include <errno.h>
|
||||
|
||||
void* allocate (int alignment, size_t bytesize)
|
||||
{
|
||||
int errorCode;
|
||||
void* ptr;
|
||||
|
||||
errorCode = posix_memalign(&ptr, alignment, bytesize);
|
||||
|
||||
if (errorCode) {
|
||||
if (errorCode == EINVAL) {
|
||||
fprintf(stderr,
|
||||
"Error: Alignment parameter is not a power of two\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
if (errorCode == ENOMEM) {
|
||||
fprintf(stderr,
|
||||
"Error: Insufficient memory to fulfill the request\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
}
|
||||
|
||||
if (ptr == NULL) {
|
||||
fprintf(stderr, "Error: posix_memalign failed!\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
return ptr;
|
||||
}
|
||||
|
||||
void* reallocate (
|
||||
void* ptr,
|
||||
int alignment,
|
||||
size_t newBytesize,
|
||||
size_t oldBytesize)
|
||||
{
|
||||
void* newarray = allocate(alignment, newBytesize);
|
||||
|
||||
if(ptr != NULL) {
|
||||
memcpy(newarray, ptr, oldBytesize);
|
||||
free(ptr);
|
||||
}
|
||||
|
||||
return newarray;
|
||||
}
|
||||
200
src/core/atom.c
Normal file
200
src/core/atom.c
Normal file
@@ -0,0 +1,200 @@
|
||||
/*
|
||||
* =======================================================================================
|
||||
*
|
||||
* Authors: Jan Eitzinger (je), jan.eitzinger@fau.de
|
||||
* Rafael Ravedutti (rr), rafaelravedutti@gmail.com
|
||||
*
|
||||
* Copyright (c) 2021 RRZE, University Erlangen-Nuremberg
|
||||
*
|
||||
* This file is part of MD-Bench.
|
||||
*
|
||||
* MD-Bench is free software: you can redistribute it and/or modify it
|
||||
* under the terms of the GNU Lesser General Public License as published
|
||||
* by the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
|
||||
* PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
|
||||
* details.
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License along
|
||||
* with MD-Bench. If not, see <https://www.gnu.org/licenses/>.
|
||||
* =======================================================================================
|
||||
*/
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
|
||||
#include <atom.h>
|
||||
#include <allocate.h>
|
||||
#include <util.h>
|
||||
|
||||
#define DELTA 20000
|
||||
|
||||
void initAtom(Atom *atom)
|
||||
{
|
||||
atom->x = NULL; atom->y = NULL; atom->z = NULL;
|
||||
atom->vx = NULL; atom->vy = NULL; atom->vz = NULL;
|
||||
atom->fx = NULL; atom->fy = NULL; atom->fz = NULL;
|
||||
atom->Natoms = 0;
|
||||
atom->Nlocal = 0;
|
||||
atom->Nghost = 0;
|
||||
atom->Nmax = 0;
|
||||
}
|
||||
|
||||
void createAtom(Atom *atom, Parameter *param)
|
||||
{
|
||||
MD_FLOAT xlo = 0.0; MD_FLOAT xhi = param->xprd;
|
||||
MD_FLOAT ylo = 0.0; MD_FLOAT yhi = param->yprd;
|
||||
MD_FLOAT zlo = 0.0; MD_FLOAT zhi = param->zprd;
|
||||
atom->Natoms = 4 * param->nx * param->ny * param->nz;
|
||||
atom->Nlocal = 0;
|
||||
MD_FLOAT alat = pow((4.0 / param->rho), (1.0 / 3.0));
|
||||
int ilo = (int) (xlo / (0.5 * alat) - 1);
|
||||
int ihi = (int) (xhi / (0.5 * alat) + 1);
|
||||
int jlo = (int) (ylo / (0.5 * alat) - 1);
|
||||
int jhi = (int) (yhi / (0.5 * alat) + 1);
|
||||
int klo = (int) (zlo / (0.5 * alat) - 1);
|
||||
int khi = (int) (zhi / (0.5 * alat) + 1);
|
||||
|
||||
ilo = MAX(ilo, 0);
|
||||
ihi = MIN(ihi, 2 * param->nx - 1);
|
||||
jlo = MAX(jlo, 0);
|
||||
jhi = MIN(jhi, 2 * param->ny - 1);
|
||||
klo = MAX(klo, 0);
|
||||
khi = MIN(khi, 2 * param->nz - 1);
|
||||
|
||||
MD_FLOAT xtmp, ytmp, ztmp, vxtmp, vytmp, vztmp;
|
||||
int i, j, k, m, n;
|
||||
int sx = 0; int sy = 0; int sz = 0;
|
||||
int ox = 0; int oy = 0; int oz = 0;
|
||||
int subboxdim = 8;
|
||||
|
||||
while(oz * subboxdim <= khi) {
|
||||
|
||||
k = oz * subboxdim + sz;
|
||||
j = oy * subboxdim + sy;
|
||||
i = ox * subboxdim + sx;
|
||||
|
||||
if(((i + j + k) % 2 == 0) &&
|
||||
(i >= ilo) && (i <= ihi) &&
|
||||
(j >= jlo) && (j <= jhi) &&
|
||||
(k >= klo) && (k <= khi)) {
|
||||
|
||||
xtmp = 0.5 * alat * i;
|
||||
ytmp = 0.5 * alat * j;
|
||||
ztmp = 0.5 * alat * k;
|
||||
|
||||
if( xtmp >= xlo && xtmp < xhi &&
|
||||
ytmp >= ylo && ytmp < yhi &&
|
||||
ztmp >= zlo && ztmp < zhi ) {
|
||||
|
||||
n = k * (2 * param->ny) * (2 * param->nx) +
|
||||
j * (2 * param->nx) +
|
||||
i + 1;
|
||||
|
||||
for(m = 0; m < 5; m++) {
|
||||
myrandom(&n);
|
||||
}
|
||||
vxtmp = myrandom(&n);
|
||||
|
||||
for(m = 0; m < 5; m++){
|
||||
myrandom(&n);
|
||||
}
|
||||
vytmp = myrandom(&n);
|
||||
|
||||
for(m = 0; m < 5; m++) {
|
||||
myrandom(&n);
|
||||
}
|
||||
vztmp = myrandom(&n);
|
||||
|
||||
if(atom->Nlocal == atom->Nmax) {
|
||||
growAtom(atom);
|
||||
}
|
||||
|
||||
atom_x(atom->Nlocal) = xtmp;
|
||||
atom_y(atom->Nlocal) = ytmp;
|
||||
atom_z(atom->Nlocal) = ztmp;
|
||||
atom->vx[atom->Nlocal] = vxtmp;
|
||||
atom->vy[atom->Nlocal] = vytmp;
|
||||
atom->vz[atom->Nlocal] = vztmp;
|
||||
atom->Nlocal++;
|
||||
}
|
||||
}
|
||||
|
||||
sx++;
|
||||
|
||||
if(sx == subboxdim) { sx = 0; sy++; }
|
||||
if(sy == subboxdim) { sy = 0; sz++; }
|
||||
if(sz == subboxdim) { sz = 0; ox++; }
|
||||
if(ox * subboxdim > ihi) { ox = 0; oy++; }
|
||||
if(oy * subboxdim > jhi) { oy = 0; oz++; }
|
||||
}
|
||||
}
|
||||
|
||||
void growAtom(Atom *atom)
|
||||
{
|
||||
int nold = atom->Nmax;
|
||||
atom->Nmax += DELTA;
|
||||
|
||||
#ifdef AOS
|
||||
atom->x = (MD_FLOAT*) reallocate(atom->x, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT) * 3, nold * sizeof(MD_FLOAT) * 3);
|
||||
#else
|
||||
atom->x = (MD_FLOAT*) reallocate(atom->x, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
||||
atom->y = (MD_FLOAT*) reallocate(atom->y, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
||||
atom->z = (MD_FLOAT*) reallocate(atom->z, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
||||
#endif
|
||||
atom->vx = (MD_FLOAT*) reallocate(atom->vx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
||||
atom->vy = (MD_FLOAT*) reallocate(atom->vy, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
||||
atom->vz = (MD_FLOAT*) reallocate(atom->vz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
||||
atom->fx = (MD_FLOAT*) reallocate(atom->fx, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
||||
atom->fy = (MD_FLOAT*) reallocate(atom->fy, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
||||
atom->fz = (MD_FLOAT*) reallocate(atom->fz, ALIGNMENT, atom->Nmax * sizeof(MD_FLOAT), nold * sizeof(MD_FLOAT));
|
||||
}
|
||||
|
||||
|
||||
/* void sortAtom() */
|
||||
/* { */
|
||||
/* binatoms(neighbor); */
|
||||
/* int* binpos = neighbor->bincount; */
|
||||
/* int* bins = neighbor->bins; */
|
||||
|
||||
/* int mbins = neighbor->mbins; */
|
||||
/* int atoms_per_bin = neighbor->atoms_per_bin; */
|
||||
|
||||
/* for(int i=1; i<mbins; i++) { */
|
||||
/* binpos[i] += binpos[i-1]; */
|
||||
/* } */
|
||||
|
||||
/* double* new_x = (double*) malloc(Nmax * sizeof(double)); */
|
||||
/* double* new_y = (double*) malloc(Nmax * sizeof(double)); */
|
||||
/* double* new_z = (double*) malloc(Nmax * sizeof(double)); */
|
||||
/* double* new_vx = (double*) malloc(Nmax * sizeof(double)); */
|
||||
/* double* new_vy = (double*) malloc(Nmax * sizeof(double)); */
|
||||
/* double* new_vz = (double*) malloc(Nmax * sizeof(double)); */
|
||||
|
||||
/* double* old_x = x; double* old_y = y; double* old_z = z; */
|
||||
/* double* old_vx = vx; double* old_vy = vy; double* old_vz = vz; */
|
||||
|
||||
/* for(int mybin = 0; mybin<mbins; mybin++) { */
|
||||
/* int start = mybin>0?binpos[mybin-1]:0; */
|
||||
/* int count = binpos[mybin] - start; */
|
||||
|
||||
/* for(int k=0; k<count; k++) { */
|
||||
/* int new_i = start + k; */
|
||||
/* int old_i = bins[mybin * atoms_per_bin + k]; */
|
||||
/* new_x[new_i] = old_x[old_i]; */
|
||||
/* new_y[new_i] = old_y[old_i]; */
|
||||
/* new_z[new_i] = old_z[old_i]; */
|
||||
/* new_vx[new_i] = old_vx[old_i]; */
|
||||
/* new_vy[new_i] = old_vy[old_i]; */
|
||||
/* new_vz[new_i] = old_vz[old_i]; */
|
||||
/* } */
|
||||
/* } */
|
||||
|
||||
/* free(x); free(y); free(z); */
|
||||
/* free(vx); free(vy); free(vz); */
|
||||
/* x = new_x; y = new_y; z = new_z; */
|
||||
/* vx = new_vx; vy = new_vy; vz = new_vz; */
|
||||
/* } */
|
||||
89
src/core/force.c
Normal file
89
src/core/force.c
Normal file
@@ -0,0 +1,89 @@
|
||||
/*
|
||||
* =======================================================================================
|
||||
*
|
||||
* Author: Jan Eitzinger (je), jan.eitzinger@fau.de
|
||||
* Copyright (c) 2021 RRZE, University Erlangen-Nuremberg
|
||||
*
|
||||
* This file is part of MD-Bench.
|
||||
*
|
||||
* MD-Bench is free software: you can redistribute it and/or modify it
|
||||
* under the terms of the GNU Lesser General Public License as published
|
||||
* by the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
|
||||
* PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
|
||||
* details.
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License along
|
||||
* with MD-Bench. If not, see <https://www.gnu.org/licenses/>.
|
||||
* =======================================================================================
|
||||
*/
|
||||
#include <likwid-marker.h>
|
||||
|
||||
#include <timing.h>
|
||||
#include <neighbor.h>
|
||||
#include <parameter.h>
|
||||
#include <atom.h>
|
||||
|
||||
double computeForce(
|
||||
Parameter *param,
|
||||
Atom *atom,
|
||||
Neighbor *neighbor)
|
||||
{
|
||||
int Nlocal = atom->Nlocal;
|
||||
int* neighs;
|
||||
MD_FLOAT cutforcesq = param->cutforce * param->cutforce;
|
||||
MD_FLOAT sigma6 = param->sigma6;
|
||||
MD_FLOAT epsilon = param->epsilon;
|
||||
MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
|
||||
MD_FLOAT S, E;
|
||||
|
||||
S = getTimeStamp();
|
||||
for(int i = 0; i < Nlocal; i++) {
|
||||
fx[i] = 0.0;
|
||||
fy[i] = 0.0;
|
||||
fz[i] = 0.0;
|
||||
}
|
||||
|
||||
LIKWID_MARKER_START("force");
|
||||
|
||||
#pragma omp parallel for
|
||||
for(int i = 0; i < Nlocal; i++) {
|
||||
neighs = &neighbor->neighbors[i * neighbor->maxneighs];
|
||||
int numneighs = neighbor->numneigh[i];
|
||||
MD_FLOAT xtmp = atom_x(i);
|
||||
MD_FLOAT ytmp = atom_y(i);
|
||||
MD_FLOAT ztmp = atom_z(i);
|
||||
|
||||
MD_FLOAT fix = 0;
|
||||
MD_FLOAT fiy = 0;
|
||||
MD_FLOAT fiz = 0;
|
||||
|
||||
for(int k = 0; k < numneighs; k++) {
|
||||
int j = neighs[k];
|
||||
MD_FLOAT delx = xtmp - atom_x(j);
|
||||
MD_FLOAT dely = ytmp - atom_y(j);
|
||||
MD_FLOAT delz = ztmp - atom_z(j);
|
||||
MD_FLOAT rsq = delx * delx + dely * dely + delz * delz;
|
||||
|
||||
if(rsq < cutforcesq) {
|
||||
MD_FLOAT sr2 = 1.0 / rsq;
|
||||
MD_FLOAT sr6 = sr2 * sr2 * sr2 * sigma6;
|
||||
MD_FLOAT force = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon;
|
||||
fix += delx * force;
|
||||
fiy += dely * force;
|
||||
fiz += delz * force;
|
||||
}
|
||||
}
|
||||
|
||||
fx[i] += fix;
|
||||
fy[i] += fiy;
|
||||
fz[i] += fiz;
|
||||
}
|
||||
LIKWID_MARKER_STOP("force");
|
||||
E = getTimeStamp();
|
||||
|
||||
return E-S;
|
||||
}
|
||||
334
src/core/neighbor.c
Normal file
334
src/core/neighbor.c
Normal file
@@ -0,0 +1,334 @@
|
||||
/*
|
||||
* =======================================================================================
|
||||
*
|
||||
* Author: Jan Eitzinger (je), jan.eitzinger@fau.de
|
||||
* Copyright (c) 2021 RRZE, University Erlangen-Nuremberg
|
||||
*
|
||||
* This file is part of MD-Bench.
|
||||
*
|
||||
* MD-Bench is free software: you can redistribute it and/or modify it
|
||||
* under the terms of the GNU Lesser General Public License as published
|
||||
* by the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
|
||||
* PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
|
||||
* details.
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License along
|
||||
* with MD-Bench. If not, see <https://www.gnu.org/licenses/>.
|
||||
* =======================================================================================
|
||||
*/
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
|
||||
#include <neighbor.h>
|
||||
#include <parameter.h>
|
||||
#include <atom.h>
|
||||
|
||||
#define SMALL 1.0e-6
|
||||
#define FACTOR 0.999
|
||||
|
||||
static MD_FLOAT xprd, yprd, zprd;
|
||||
static MD_FLOAT bininvx, bininvy, bininvz;
|
||||
static int mbinxlo, mbinylo, mbinzlo;
|
||||
static int nbinx, nbiny, nbinz;
|
||||
static int mbinx, mbiny, mbinz; // n bins in x, y, z
|
||||
static int *bincount;
|
||||
static int *bins;
|
||||
static int mbins; //total number of bins
|
||||
static int atoms_per_bin; // max atoms per bin
|
||||
static MD_FLOAT cutneigh;
|
||||
static MD_FLOAT cutneighsq; // neighbor cutoff squared
|
||||
static int nmax;
|
||||
static int nstencil; // # of bins in stencil
|
||||
static int* stencil; // stencil list of bin offsets
|
||||
static MD_FLOAT binsizex, binsizey, binsizez;
|
||||
|
||||
static int coord2bin(MD_FLOAT, MD_FLOAT , MD_FLOAT);
|
||||
static MD_FLOAT bindist(int, int, int);
|
||||
|
||||
/* exported subroutines */
|
||||
void initNeighbor(Neighbor *neighbor, Parameter *param)
|
||||
{
|
||||
MD_FLOAT neighscale = 5.0 / 6.0;
|
||||
xprd = param->nx * param->lattice;
|
||||
yprd = param->ny * param->lattice;
|
||||
zprd = param->nz * param->lattice;
|
||||
cutneigh = param->cutneigh;
|
||||
nbinx = neighscale * param->nx;
|
||||
nbiny = neighscale * param->ny;
|
||||
nbinz = neighscale * param->nz;
|
||||
nmax = 0;
|
||||
atoms_per_bin = 8;
|
||||
stencil = NULL;
|
||||
bins = NULL;
|
||||
bincount = NULL;
|
||||
neighbor->maxneighs = 100;
|
||||
neighbor->numneigh = NULL;
|
||||
neighbor->neighbors = NULL;
|
||||
}
|
||||
|
||||
void setupNeighbor()
|
||||
{
|
||||
MD_FLOAT coord;
|
||||
int mbinxhi, mbinyhi, mbinzhi;
|
||||
int nextx, nexty, nextz;
|
||||
MD_FLOAT xlo = 0.0; MD_FLOAT xhi = xprd;
|
||||
MD_FLOAT ylo = 0.0; MD_FLOAT yhi = yprd;
|
||||
MD_FLOAT zlo = 0.0; MD_FLOAT zhi = zprd;
|
||||
|
||||
cutneighsq = cutneigh * cutneigh;
|
||||
binsizex = xprd / nbinx;
|
||||
binsizey = yprd / nbiny;
|
||||
binsizez = zprd / nbinz;
|
||||
bininvx = 1.0 / binsizex;
|
||||
bininvy = 1.0 / binsizey;
|
||||
bininvz = 1.0 / binsizez;
|
||||
|
||||
coord = xlo - cutneigh - SMALL * xprd;
|
||||
mbinxlo = (int) (coord * bininvx);
|
||||
if (coord < 0.0) {
|
||||
mbinxlo = mbinxlo - 1;
|
||||
}
|
||||
coord = xhi + cutneigh + SMALL * xprd;
|
||||
mbinxhi = (int) (coord * bininvx);
|
||||
|
||||
coord = ylo - cutneigh - SMALL * yprd;
|
||||
mbinylo = (int) (coord * bininvy);
|
||||
if (coord < 0.0) {
|
||||
mbinylo = mbinylo - 1;
|
||||
}
|
||||
coord = yhi + cutneigh + SMALL * yprd;
|
||||
mbinyhi = (int) (coord * bininvy);
|
||||
|
||||
coord = zlo - cutneigh - SMALL * zprd;
|
||||
mbinzlo = (int) (coord * bininvz);
|
||||
if (coord < 0.0) {
|
||||
mbinzlo = mbinzlo - 1;
|
||||
}
|
||||
coord = zhi + cutneigh + SMALL * zprd;
|
||||
mbinzhi = (int) (coord * bininvz);
|
||||
|
||||
mbinxlo = mbinxlo - 1;
|
||||
mbinxhi = mbinxhi + 1;
|
||||
mbinx = mbinxhi - mbinxlo + 1;
|
||||
|
||||
mbinylo = mbinylo - 1;
|
||||
mbinyhi = mbinyhi + 1;
|
||||
mbiny = mbinyhi - mbinylo + 1;
|
||||
|
||||
mbinzlo = mbinzlo - 1;
|
||||
mbinzhi = mbinzhi + 1;
|
||||
mbinz = mbinzhi - mbinzlo + 1;
|
||||
|
||||
nextx = (int) (cutneigh * bininvx);
|
||||
if(nextx * binsizex < FACTOR * cutneigh) nextx++;
|
||||
|
||||
nexty = (int) (cutneigh * bininvy);
|
||||
if(nexty * binsizey < FACTOR * cutneigh) nexty++;
|
||||
|
||||
nextz = (int) (cutneigh * bininvz);
|
||||
if(nextz * binsizez < FACTOR * cutneigh) nextz++;
|
||||
|
||||
if (stencil) {
|
||||
free(stencil);
|
||||
}
|
||||
|
||||
stencil = (int*) malloc(
|
||||
(2 * nextz + 1) * (2 * nexty + 1) * (2 * nextx + 1) * sizeof(int));
|
||||
|
||||
nstencil = 0;
|
||||
int kstart = -nextz;
|
||||
|
||||
for(int k = kstart; k <= nextz; k++) {
|
||||
for(int j = -nexty; j <= nexty; j++) {
|
||||
for(int i = -nextx; i <= nextx; i++) {
|
||||
if(bindist(i, j, k) < cutneighsq) {
|
||||
stencil[nstencil++] =
|
||||
k * mbiny * mbinx + j * mbinx + i;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
mbins = mbinx * mbiny * mbinz;
|
||||
|
||||
if (bincount) {
|
||||
free(bincount);
|
||||
}
|
||||
bincount = (int*) malloc(mbins * sizeof(int));
|
||||
|
||||
if (bins) {
|
||||
free(bins);
|
||||
}
|
||||
bins = (int*) malloc(mbins * atoms_per_bin * sizeof(int));
|
||||
}
|
||||
|
||||
void buildNeighbor(Atom *atom, Neighbor *neighbor)
|
||||
{
|
||||
int nall = atom->Nlocal + atom->Nghost;
|
||||
|
||||
/* extend atom arrays if necessary */
|
||||
if(nall > nmax) {
|
||||
nmax = nall;
|
||||
if(neighbor->numneigh) free(neighbor->numneigh);
|
||||
if(neighbor->neighbors) free(neighbor->neighbors);
|
||||
neighbor->numneigh = (int*) malloc(nmax * sizeof(int));
|
||||
neighbor->neighbors = (int*) malloc(nmax * neighbor->maxneighs * sizeof(int*));
|
||||
}
|
||||
|
||||
/* bin local & ghost atoms */
|
||||
binatoms(atom);
|
||||
int resize = 1;
|
||||
|
||||
/* loop over each atom, storing neighbors */
|
||||
while(resize) {
|
||||
int new_maxneighs = neighbor->maxneighs;
|
||||
resize = 0;
|
||||
|
||||
for(int i = 0; i < atom->Nlocal; i++) {
|
||||
int* neighptr = &(neighbor->neighbors[i * neighbor->maxneighs]);
|
||||
int n = 0;
|
||||
MD_FLOAT xtmp = atom_x(i);
|
||||
MD_FLOAT ytmp = atom_y(i);
|
||||
MD_FLOAT ztmp = atom_z(i);
|
||||
int ibin = coord2bin(xtmp, ytmp, ztmp);
|
||||
|
||||
for(int k = 0; k < nstencil; k++) {
|
||||
int jbin = ibin + stencil[k];
|
||||
int* loc_bin = &bins[jbin * atoms_per_bin];
|
||||
|
||||
for(int m = 0; m < bincount[jbin]; m++) {
|
||||
int j = loc_bin[m];
|
||||
|
||||
if ( j == i ){
|
||||
continue;
|
||||
}
|
||||
|
||||
MD_FLOAT delx = xtmp - atom_x(j);
|
||||
MD_FLOAT dely = ytmp - atom_y(j);
|
||||
MD_FLOAT delz = ztmp - atom_z(j);
|
||||
MD_FLOAT rsq = delx * delx + dely * dely + delz * delz;
|
||||
|
||||
if( rsq <= cutneighsq ) {
|
||||
neighptr[n++] = j;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
neighbor->numneigh[i] = n;
|
||||
|
||||
if(n >= neighbor->maxneighs) {
|
||||
resize = 1;
|
||||
|
||||
if(n >= new_maxneighs) {
|
||||
new_maxneighs = n;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if(resize) {
|
||||
printf("RESIZE %d\n", neighbor->maxneighs);
|
||||
neighbor->maxneighs = new_maxneighs * 1.2;
|
||||
free(neighbor->neighbors);
|
||||
neighbor->neighbors = (int*) malloc(atom->Nmax * neighbor->maxneighs * sizeof(int));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* internal subroutines */
|
||||
MD_FLOAT bindist(int i, int j, int k)
|
||||
{
|
||||
MD_FLOAT delx, dely, delz;
|
||||
|
||||
if(i > 0) {
|
||||
delx = (i - 1) * binsizex;
|
||||
} else if(i == 0) {
|
||||
delx = 0.0;
|
||||
} else {
|
||||
delx = (i + 1) * binsizex;
|
||||
}
|
||||
|
||||
if(j > 0) {
|
||||
dely = (j - 1) * binsizey;
|
||||
} else if(j == 0) {
|
||||
dely = 0.0;
|
||||
} else {
|
||||
dely = (j + 1) * binsizey;
|
||||
}
|
||||
|
||||
if(k > 0) {
|
||||
delz = (k - 1) * binsizez;
|
||||
} else if(k == 0) {
|
||||
delz = 0.0;
|
||||
} else {
|
||||
delz = (k + 1) * binsizez;
|
||||
}
|
||||
|
||||
return (delx * delx + dely * dely + delz * delz);
|
||||
}
|
||||
|
||||
int coord2bin(MD_FLOAT xin, MD_FLOAT yin, MD_FLOAT zin)
|
||||
{
|
||||
int ix, iy, iz;
|
||||
|
||||
if(xin >= xprd) {
|
||||
ix = (int)((xin - xprd) * bininvx) + nbinx - mbinxlo;
|
||||
} else if(xin >= 0.0) {
|
||||
ix = (int)(xin * bininvx) - mbinxlo;
|
||||
} else {
|
||||
ix = (int)(xin * bininvx) - mbinxlo - 1;
|
||||
}
|
||||
|
||||
if(yin >= yprd) {
|
||||
iy = (int)((yin - yprd) * bininvy) + nbiny - mbinylo;
|
||||
} else if(yin >= 0.0) {
|
||||
iy = (int)(yin * bininvy) - mbinylo;
|
||||
} else {
|
||||
iy = (int)(yin * bininvy) - mbinylo - 1;
|
||||
}
|
||||
|
||||
if(zin >= zprd) {
|
||||
iz = (int)((zin - zprd) * bininvz) + nbinz - mbinzlo;
|
||||
} else if(zin >= 0.0) {
|
||||
iz = (int)(zin * bininvz) - mbinzlo;
|
||||
} else {
|
||||
iz = (int)(zin * bininvz) - mbinzlo - 1;
|
||||
}
|
||||
|
||||
return (iz * mbiny * mbinx + iy * mbinx + ix + 1);
|
||||
}
|
||||
|
||||
void binatoms(Atom *atom)
|
||||
{
|
||||
int nall = atom->Nlocal + atom->Nghost;
|
||||
int resize = 1;
|
||||
|
||||
while(resize > 0) {
|
||||
resize = 0;
|
||||
|
||||
for(int i = 0; i < mbins; i++) {
|
||||
bincount[i] = 0;
|
||||
}
|
||||
|
||||
for(int i = 0; i < nall; i++) {
|
||||
int ibin = coord2bin(atom_x(i), atom_y(i), atom_z(i));
|
||||
|
||||
if(bincount[ibin] < atoms_per_bin) {
|
||||
int ac = bincount[ibin]++;
|
||||
bins[ibin * atoms_per_bin + ac] = i;
|
||||
} else {
|
||||
resize = 1;
|
||||
}
|
||||
}
|
||||
|
||||
if(resize) {
|
||||
free(bins);
|
||||
atoms_per_bin *= 2;
|
||||
bins = (int*) malloc(mbins * atoms_per_bin * sizeof(int));
|
||||
}
|
||||
}
|
||||
}
|
||||
168
src/core/pbc.c
Normal file
168
src/core/pbc.c
Normal file
@@ -0,0 +1,168 @@
|
||||
/*
|
||||
* =======================================================================================
|
||||
*
|
||||
* Author: Jan Eitzinger (je), jan.eitzinger@fau.de
|
||||
* Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
|
||||
*
|
||||
* This file is part of MD-Bench.
|
||||
*
|
||||
* MD-Bench is free software: you can redistribute it and/or modify it
|
||||
* under the terms of the GNU Lesser General Public License as published
|
||||
* by the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
|
||||
* PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
|
||||
* details.
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License along
|
||||
* with MD-Bench. If not, see <https://www.gnu.org/licenses/>.
|
||||
* =======================================================================================
|
||||
*/
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
|
||||
#include <pbc.h>
|
||||
#include <atom.h>
|
||||
#include <allocate.h>
|
||||
|
||||
#define DELTA 20000
|
||||
|
||||
static int NmaxGhost;
|
||||
static int *BorderMap;
|
||||
static int *PBCx, *PBCy, *PBCz;
|
||||
|
||||
static void growPbc();
|
||||
|
||||
/* exported subroutines */
|
||||
void initPbc()
|
||||
{
|
||||
NmaxGhost = 0;
|
||||
BorderMap = NULL;
|
||||
PBCx = NULL; PBCy = NULL; PBCz = NULL;
|
||||
}
|
||||
|
||||
/* update coordinates of ghost atoms */
|
||||
/* uses mapping created in setupPbc */
|
||||
void updatePbc(Atom *atom, Parameter *param)
|
||||
{
|
||||
int nlocal = atom->Nlocal;
|
||||
MD_FLOAT xprd = param->xprd;
|
||||
MD_FLOAT yprd = param->yprd;
|
||||
MD_FLOAT zprd = param->zprd;
|
||||
|
||||
for(int i = 0; i < atom->Nghost; i++) {
|
||||
atom_x(nlocal + i) = atom_x(BorderMap[i]) + PBCx[i] * xprd;
|
||||
atom_y(nlocal + i) = atom_y(BorderMap[i]) + PBCy[i] * yprd;
|
||||
atom_z(nlocal + i) = atom_z(BorderMap[i]) + PBCz[i] * zprd;
|
||||
}
|
||||
}
|
||||
|
||||
/* relocate atoms that have left domain according
|
||||
* to periodic boundary conditions */
|
||||
void updateAtomsPbc(Atom *atom, Parameter *param)
|
||||
{
|
||||
MD_FLOAT xprd = param->xprd;
|
||||
MD_FLOAT yprd = param->yprd;
|
||||
MD_FLOAT zprd = param->zprd;
|
||||
|
||||
for(int i = 0; i < atom->Nlocal; i++) {
|
||||
|
||||
if(atom_x(i) < 0.0) {
|
||||
atom_x(i) += xprd;
|
||||
} else if(atom_x(i) >= xprd) {
|
||||
atom_x(i) -= xprd;
|
||||
}
|
||||
|
||||
if(atom_y(i) < 0.0) {
|
||||
atom_y(i) += yprd;
|
||||
} else if(atom_y(i) >= yprd) {
|
||||
atom_y(i) -= yprd;
|
||||
}
|
||||
|
||||
if(atom_z(i) < 0.0) {
|
||||
atom_z(i) += zprd;
|
||||
} else if(atom_z(i) >= zprd) {
|
||||
atom_z(i) -= zprd;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* setup periodic boundary conditions by
|
||||
* defining ghost atoms around domain
|
||||
* only creates mapping and coordinate corrections
|
||||
* that are then enforced in updatePbc */
|
||||
#define ADDGHOST(dx,dy,dz) \
|
||||
Nghost++; \
|
||||
BorderMap[Nghost] = i; \
|
||||
PBCx[Nghost] = dx; \
|
||||
PBCy[Nghost] = dy; \
|
||||
PBCz[Nghost] = dz;
|
||||
void setupPbc(Atom *atom, Parameter *param)
|
||||
{
|
||||
MD_FLOAT xprd = param->xprd;
|
||||
MD_FLOAT yprd = param->yprd;
|
||||
MD_FLOAT zprd = param->zprd;
|
||||
MD_FLOAT Cutneigh = param->cutneigh;
|
||||
int Nghost = -1;
|
||||
|
||||
for(int i = 0; i < atom->Nlocal; i++) {
|
||||
|
||||
if (atom->Nlocal + Nghost + 7 >= atom->Nmax) {
|
||||
growAtom(atom);
|
||||
}
|
||||
if (Nghost + 7 >= NmaxGhost) {
|
||||
growPbc();
|
||||
}
|
||||
|
||||
MD_FLOAT x = atom_x(i);
|
||||
MD_FLOAT y = atom_y(i);
|
||||
MD_FLOAT z = atom_z(i);
|
||||
|
||||
/* Setup ghost atoms */
|
||||
/* 6 planes */
|
||||
if (x < Cutneigh) { ADDGHOST(+1,0,0); }
|
||||
if (x >= (xprd-Cutneigh)) { ADDGHOST(-1,0,0); }
|
||||
if (y < Cutneigh) { ADDGHOST(0,+1,0); }
|
||||
if (y >= (yprd-Cutneigh)) { ADDGHOST(0,-1,0); }
|
||||
if (z < Cutneigh) { ADDGHOST(0,0,+1); }
|
||||
if (z >= (zprd-Cutneigh)) { ADDGHOST(0,0,-1); }
|
||||
/* 8 corners */
|
||||
if (x < Cutneigh && y < Cutneigh && z < Cutneigh) { ADDGHOST(+1,+1,+1); }
|
||||
if (x < Cutneigh && y >= (yprd-Cutneigh) && z < Cutneigh) { ADDGHOST(+1,-1,+1); }
|
||||
if (x < Cutneigh && y >= Cutneigh && z >= (zprd-Cutneigh)) { ADDGHOST(+1,+1,-1); }
|
||||
if (x < Cutneigh && y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(+1,-1,-1); }
|
||||
if (x >= (xprd-Cutneigh) && y < Cutneigh && z < Cutneigh) { ADDGHOST(-1,+1,+1); }
|
||||
if (x >= (xprd-Cutneigh) && y >= (yprd-Cutneigh) && z < Cutneigh) { ADDGHOST(-1,-1,+1); }
|
||||
if (x >= (xprd-Cutneigh) && y < Cutneigh && z >= (zprd-Cutneigh)) { ADDGHOST(-1,+1,-1); }
|
||||
if (x >= (xprd-Cutneigh) && y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(-1,-1,-1); }
|
||||
/* 12 edges */
|
||||
if (x < Cutneigh && z < Cutneigh) { ADDGHOST(+1,0,+1); }
|
||||
if (x < Cutneigh && z >= (zprd-Cutneigh)) { ADDGHOST(+1,0,-1); }
|
||||
if (x >= (xprd-Cutneigh) && z < Cutneigh) { ADDGHOST(-1,0,+1); }
|
||||
if (x >= (xprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(-1,0,-1); }
|
||||
if (y < Cutneigh && z < Cutneigh) { ADDGHOST(0,+1,+1); }
|
||||
if (y < Cutneigh && z >= (zprd-Cutneigh)) { ADDGHOST(0,+1,-1); }
|
||||
if (y >= (yprd-Cutneigh) && z < Cutneigh) { ADDGHOST(0,-1,+1); }
|
||||
if (y >= (yprd-Cutneigh) && z >= (zprd-Cutneigh)) { ADDGHOST(0,-1,-1); }
|
||||
if (y < Cutneigh && x < Cutneigh) { ADDGHOST(+1,+1,0); }
|
||||
if (y < Cutneigh && x >= (xprd-Cutneigh)) { ADDGHOST(-1,+1,0); }
|
||||
if (y >= (yprd-Cutneigh) && x < Cutneigh) { ADDGHOST(+1,-1,0); }
|
||||
if (y >= (yprd-Cutneigh) && x >= (xprd-Cutneigh)) { ADDGHOST(-1,-1,0); }
|
||||
}
|
||||
// increase by one to make it the ghost atom count
|
||||
atom->Nghost = Nghost + 1;
|
||||
}
|
||||
|
||||
/* internal subroutines */
|
||||
void growPbc()
|
||||
{
|
||||
int nold = NmaxGhost;
|
||||
NmaxGhost += DELTA;
|
||||
|
||||
BorderMap = (int*) reallocate(BorderMap, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));
|
||||
PBCx = (int*) reallocate(PBCx, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));
|
||||
PBCy = (int*) reallocate(PBCy, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));
|
||||
PBCz = (int*) reallocate(PBCz, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int));
|
||||
}
|
||||
129
src/core/thermo.c
Normal file
129
src/core/thermo.c
Normal file
@@ -0,0 +1,129 @@
|
||||
/*
|
||||
* =======================================================================================
|
||||
*
|
||||
* Author: Jan Eitzinger (je), jan.eitzinger@fau.de
|
||||
* Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
|
||||
*
|
||||
* This file is part of MD-Bench.
|
||||
*
|
||||
* MD-Bench is free software: you can redistribute it and/or modify it
|
||||
* under the terms of the GNU Lesser General Public License as published
|
||||
* by the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
|
||||
* PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
|
||||
* details.
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License along
|
||||
* with MD-Bench. If not, see <https://www.gnu.org/licenses/>.
|
||||
* =======================================================================================
|
||||
*/
|
||||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
|
||||
#include <thermo.h>
|
||||
|
||||
static int *steparr;
|
||||
static MD_FLOAT *tmparr;
|
||||
static MD_FLOAT *engarr;
|
||||
static MD_FLOAT *prsarr;
|
||||
static MD_FLOAT mvv2e;
|
||||
static int dof_boltz;
|
||||
static MD_FLOAT t_scale;
|
||||
static MD_FLOAT p_scale;
|
||||
static MD_FLOAT e_scale;
|
||||
static MD_FLOAT t_act;
|
||||
static MD_FLOAT p_act;
|
||||
static MD_FLOAT e_act;
|
||||
static int mstat;
|
||||
|
||||
/* exported subroutines */
|
||||
void setupThermo(Parameter *param, int natoms)
|
||||
{
|
||||
int maxstat = param->ntimes / param->nstat + 2;
|
||||
|
||||
steparr = (int*) malloc(maxstat * sizeof(int));
|
||||
tmparr = (MD_FLOAT*) malloc(maxstat * sizeof(MD_FLOAT));
|
||||
engarr = (MD_FLOAT*) malloc(maxstat * sizeof(MD_FLOAT));
|
||||
prsarr = (MD_FLOAT*) malloc(maxstat * sizeof(MD_FLOAT));
|
||||
|
||||
mvv2e = 1.0;
|
||||
dof_boltz = (natoms * 3 - 3);
|
||||
t_scale = mvv2e / dof_boltz;
|
||||
p_scale = 1.0 / 3 / param->xprd / param->yprd / param->zprd;
|
||||
e_scale = 0.5;
|
||||
|
||||
printf("step\ttemp\t\tpressure\n");
|
||||
}
|
||||
|
||||
void computeThermo(int iflag, Parameter *param, Atom *atom)
|
||||
{
|
||||
MD_FLOAT t = 0.0, p;
|
||||
MD_FLOAT* vx = atom->vx;
|
||||
MD_FLOAT* vy = atom->vy;
|
||||
MD_FLOAT* vz = atom->vz;
|
||||
|
||||
for(int i = 0; i < atom->Nlocal; i++) {
|
||||
t += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) * param->mass;
|
||||
}
|
||||
|
||||
t = t * t_scale;
|
||||
p = (t * dof_boltz) * p_scale;
|
||||
|
||||
int istep = iflag;
|
||||
|
||||
if(iflag == -1){
|
||||
istep = param->ntimes;
|
||||
}
|
||||
if(iflag == 0){
|
||||
mstat = 0;
|
||||
}
|
||||
|
||||
steparr[mstat] = istep;
|
||||
tmparr[mstat] = t;
|
||||
prsarr[mstat] = p;
|
||||
mstat++;
|
||||
fprintf(stdout, "%i\t%e\t%e\n", istep, t, p);
|
||||
}
|
||||
|
||||
void adjustThermo(Parameter *param, Atom *atom)
|
||||
{
|
||||
/* zero center-of-mass motion */
|
||||
MD_FLOAT vxtot = 0.0; MD_FLOAT vytot = 0.0; MD_FLOAT vztot = 0.0;
|
||||
MD_FLOAT* vx = atom->vx; MD_FLOAT* vy = atom->vy; MD_FLOAT* vz = atom->vz;
|
||||
|
||||
for(int i = 0; i < atom->Nlocal; i++) {
|
||||
vxtot += vx[i];
|
||||
vytot += vy[i];
|
||||
vztot += vz[i];
|
||||
}
|
||||
|
||||
vxtot = vxtot / atom->Natoms;
|
||||
vytot = vytot / atom->Natoms;
|
||||
vztot = vztot / atom->Natoms;
|
||||
|
||||
for(int i = 0; i < atom->Nlocal; i++) {
|
||||
vx[i] -= vxtot;
|
||||
vy[i] -= vytot;
|
||||
vz[i] -= vztot;
|
||||
}
|
||||
|
||||
t_act = 0;
|
||||
MD_FLOAT t = 0.0;
|
||||
|
||||
for(int i = 0; i < atom->Nlocal; i++) {
|
||||
t += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) * param->mass;
|
||||
}
|
||||
|
||||
t *= t_scale;
|
||||
MD_FLOAT factor = sqrt(param->temp / t);
|
||||
|
||||
for(int i = 0; i < atom->Nlocal; i++) {
|
||||
vx[i] *= factor;
|
||||
vy[i] *= factor;
|
||||
vz[i] *= factor;
|
||||
}
|
||||
}
|
||||
43
src/core/timing.c
Normal file
43
src/core/timing.c
Normal file
@@ -0,0 +1,43 @@
|
||||
/*
|
||||
* =======================================================================================
|
||||
*
|
||||
* Author: Jan Eitzinger (je), jan.eitzinger@fau.de
|
||||
* Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
|
||||
*
|
||||
* This file is part of MD-Bench.
|
||||
*
|
||||
* MD-Bench is free software: you can redistribute it and/or modify it
|
||||
* under the terms of the GNU Lesser General Public License as published
|
||||
* by the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
|
||||
* PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
|
||||
* details.
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License along
|
||||
* with MD-Bench. If not, see <https://www.gnu.org/licenses/>.
|
||||
* =======================================================================================
|
||||
*/
|
||||
#include <stdlib.h>
|
||||
#include <time.h>
|
||||
|
||||
double getTimeStamp()
|
||||
{
|
||||
struct timespec ts;
|
||||
clock_gettime(CLOCK_MONOTONIC, &ts);
|
||||
return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
|
||||
}
|
||||
|
||||
double getTimeResolution()
|
||||
{
|
||||
struct timespec ts;
|
||||
clock_getres(CLOCK_MONOTONIC, &ts);
|
||||
return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
|
||||
}
|
||||
|
||||
double getTimeStamp_()
|
||||
{
|
||||
return getTimeStamp();
|
||||
}
|
||||
42
src/core/util.c
Normal file
42
src/core/util.c
Normal file
@@ -0,0 +1,42 @@
|
||||
/*
|
||||
* =======================================================================================
|
||||
*
|
||||
* Author: Jan Eitzinger (je), jan.eitzinger@fau.de
|
||||
* Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
|
||||
*
|
||||
* This file is part of MD-Bench.
|
||||
*
|
||||
* MD-Bench is free software: you can redistribute it and/or modify it
|
||||
* under the terms of the GNU Lesser General Public License as published
|
||||
* by the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
|
||||
* PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
|
||||
* details.
|
||||
*
|
||||
* You should have received a copy of the GNU Lesser General Public License along
|
||||
* with MD-Bench. If not, see <https://www.gnu.org/licenses/>.
|
||||
* =======================================================================================
|
||||
*/
|
||||
#include <util.h>
|
||||
|
||||
/* Park/Miller RNG w/out MASKING, so as to be like f90s version */
|
||||
#define IA 16807
|
||||
#define IM 2147483647
|
||||
#define AM (1.0/IM)
|
||||
#define IQ 127773
|
||||
#define IR 2836
|
||||
#define MASK 123459876
|
||||
|
||||
double myrandom(int* idum)
|
||||
{
|
||||
int k= (*idum) / IQ;
|
||||
double ans;
|
||||
|
||||
*idum = IA * (*idum - k * IQ) - IR * k;
|
||||
if(*idum < 0) *idum += IM;
|
||||
ans = AM * (*idum);
|
||||
return ans;
|
||||
}
|
||||
Reference in New Issue
Block a user