Page MenuHomec4science
Files F105496066

setup.cpp
No OneTemporary
Actions

File Metadata

Created
Mon, Mar 17, 16:13

setup.cpp
View Options

/*
//@HEADER
// ************************************************************************
//
// Kokkos v. 2.0
// Copyright (2014) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact H. Carter Edwards (hcedwar@sandia.gov)
//
// ************************************************************************
//@HEADER
*/
#include <system.h>
#include <cmath>
#include <cstdio>
#include <cstdlib>
/* initialize atoms on fcc lattice in parallel fashion */
#define MAX(a,b) (a>b?a:b)
#define MIN(a,b) (a<b?a:b)
int create_system(System &system, int nx, int ny, int nz, double rho)
{
/* Box Setup */
double lattice = pow((4.0 / rho), (1.0 / 3.0));
system.box.xprd = nx * lattice;
system.box.yprd = ny * lattice;
system.box.zprd = nz * lattice;
system.box.xlo = 0;
system.box.ylo = 0;
system.box.zlo = 0;
system.box.xhi = system.box.xprd;
system.box.yhi = system.box.yprd;
system.box.zhi = system.box.zprd;
int ghost_dist = int(system.neigh_cut/lattice) + 1;
/* total # of atoms */
system.nlocal = 4 * nx * ny * nz;
system.nghost = 4 * (nx + 2 * ghost_dist) *
(ny + 2 * ghost_dist) *
(nz + 2 * ghost_dist) -
system.nlocal;
system.natoms = system.nlocal + system.nghost;
system.d_x = t_x_array("X",system.natoms);
system.h_x = Kokkos::create_mirror_view(system.d_x);
system.f = t_f_array("F",system.natoms);
/* determine loop bounds of lattice subsection that overlaps my sub-box
insure loop bounds do not exceed nx,ny,nz */
double alat = pow((4.0 / rho), (1.0 / 3.0));
int ilo = static_cast<int>(system.box.xlo / (0.5 * alat) - 1);
int ihi = static_cast<int>(system.box.xhi / (0.5 * alat) + 1);
int jlo = static_cast<int>(system.box.ylo / (0.5 * alat) - 1);
int jhi = static_cast<int>(system.box.yhi / (0.5 * alat) + 1);
int klo = static_cast<int>(system.box.zlo / (0.5 * alat) - 1);
int khi = static_cast<int>(system.box.zhi / (0.5 * alat) + 1);
ilo = MAX(ilo, 0);
ihi = MIN(ihi, 2 * nx - 1);
jlo = MAX(jlo, 0);
jhi = MIN(jhi, 2 * ny - 1);
klo = MAX(klo, 0);
khi = MIN(khi, 2 * nz - 1);
/* generates positions of atoms on fcc sublattice*/
srand(3718273);
/* create non-ghost atoms */
{
double xtmp, ytmp, ztmp;
int sx = 0;
int sy = 0;
int sz = 0;
int ox = 0;
int oy = 0;
int oz = 0;
int subboxdim = 8;
int n = 0;
int iflag = 0;
while(oz * subboxdim <= khi) {
const int k = oz * subboxdim + sz;
const int j = oy * subboxdim + sy;
const int i = ox * subboxdim + sx;
if(iflag) continue;
if(((i + j + k) % 2 == 0) &&
(i >= ilo) && (i <= ihi) &&
(j >= jlo) && (j <= jhi) &&
(k >= klo) && (k <= khi)) {
const int nold = n;
while(nold == n) {
xtmp = 0.5 * alat * i + system.delta/1000*(rand()%1000-500);
ytmp = 0.5 * alat * j + system.delta/1000*(rand()%1000-500);
ztmp = 0.5 * alat * k + system.delta/1000*(rand()%1000-500);
if(xtmp >= system.box.xlo && xtmp < system.box.xhi &&
ytmp >= system.box.ylo && ytmp < system.box.yhi &&
ztmp >= system.box.zlo && ztmp < system.box.zhi) {
system.h_x(n,0) = xtmp;
system.h_x(n,1) = ytmp;
system.h_x(n,2) = ztmp;
n++;
}
}
}
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++;
}
}
/* check that correct # of atoms were created */
if(system.nlocal != n) {
printf("Created incorrect # of atoms\n");
return 1;
}
}
/* create ghost atoms */
{
double xtmp, ytmp, ztmp;
int ilo_g = ilo - 2 * ghost_dist;
int jlo_g = jlo - 2 * ghost_dist;
int klo_g = klo - 2 * ghost_dist;
int ihi_g = ihi + 2 * ghost_dist;
int jhi_g = jhi + 2 * ghost_dist;
int khi_g = khi + 2 * ghost_dist;
int subboxdim = 8;
int sx = 0;
int sy = 0;
int sz = 0;
int ox = subboxdim * ilo_g;
int oy = subboxdim * jlo_g;
int oz = subboxdim * klo_g;
int n = system.nlocal;
int iflag = 0;
while(oz * subboxdim <= khi_g) {
const int k = oz * subboxdim + sz;
const int j = oy * subboxdim + sy;
const int i = ox * subboxdim + sx;
if(iflag) continue;
if(((i + j + k) % 2 == 0) &&
(i >= ilo_g) && (i <= ihi_g) &&
(j >= jlo_g) && (j <= jhi_g) &&
(k >= klo_g) && (k <= khi_g) &&
((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;
system.h_x(n,0) = xtmp + system.delta/1000*(rand()%1000-500);;
system.h_x(n,1) = ytmp + system.delta/1000*(rand()%1000-500);;
system.h_x(n,2) = ztmp + system.delta/1000*(rand()%1000-500);;
n++;
}
sx++;
if(sx == subboxdim) {
sx = 0;
sy++;
}
if(sy == subboxdim) {
sy = 0;
sz++;
}
if(sz == subboxdim) {
sz = 0;
ox++;
//printf("%i %i %i // %i %i %i\n",ox,oy,oz,i,j,k);
}
if(ox * subboxdim > ihi_g) {
ox = subboxdim * ilo_g;
oy++;
}
if(oy * subboxdim > jhi_g) {
oy = subboxdim * jlo_g;
oz++;
}
}
}
Kokkos::deep_copy(system.d_x,system.h_x);
return 0;
}

Event Timeline

c4science · Help