From b6b9596aecc346f53ac5fa188102a8780b723655 Mon Sep 17 00:00:00 2001
From: Rafael Ravedutti <rafaelravedutti@gmail.com>
Date: Tue, 28 Nov 2023 01:22:53 +0100
Subject: [PATCH] Setup copper lattice without reading from file
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
---
examples/lj.py | 11 +-
runtime/copper_fcc_lattice.hpp | 154 ++++++++++++++++++++++++++++
runtime/thermo.hpp | 67 ++++++++++--
src/pairs/code_gen/cgen.py | 1 +
src/pairs/sim/copper_fcc_lattice.py | 34 ++++++
src/pairs/sim/simulation.py | 16 ++-
src/pairs/sim/thermo.py | 19 ++++
7 files changed, 291 insertions(+), 11 deletions(-)
create mode 100644 runtime/copper_fcc_lattice.hpp
create mode 100644 src/pairs/sim/copper_fcc_lattice.py
create mode 100644 src/pairs/sim/thermo.py
diff --git a/examples/lj.py b/examples/lj.py
index b5986e5..4085490 100644
--- a/examples/lj.py
+++ b/examples/lj.py
@@ -29,6 +29,11 @@ ntypes = 4
sigma = 1.0
epsilon = 1.0
sigma6 = sigma ** 6
+nx = 32
+ny = 32
+nz = 32
+rho = 0.8442
+temp = 1.44
psim = pairs.simulation("lj", [pairs.point_mass()], timesteps=200, double_prec=True)
psim.add_position('position')
@@ -38,12 +43,12 @@ psim.add_property('force', pairs.vector(), volatile=True)
psim.add_feature('type', ntypes)
psim.add_feature_property('type', 'epsilon', pairs.real(), [sigma for i in range(ntypes * ntypes)])
psim.add_feature_property('type', 'sigma6', pairs.real(), [epsilon for i in range(ntypes * ntypes)])
-psim.set_domain([0.0, 0.0, 0.0, 53.747078, 53.747078, 53.747078])
-psim.read_particle_data("data/minimd_setup_32x32x32.input", ['type', 'mass', 'position', 'linear_velocity', 'flags'], pairs.point_mass())
+psim.copper_fcc_lattice(nx, ny, nz, rho, temp, ntypes)
psim.reneighbor_every(20)
#psim.compute_half()
psim.build_neighbor_lists(cutoff_radius + skin)
-psim.vtk_output(f"output/lj_{target}")
+psim.compute_thermo(100)
+#psim.vtk_output(f"output/lj_{target}")
psim.compute(initial_integrate, symbols={'dt': dt}, pre_step=True, skip_first=True)
psim.compute(lj, cutoff_radius)
psim.compute(final_integrate, symbols={'dt': dt}, skip_first=True)
diff --git a/runtime/copper_fcc_lattice.hpp b/runtime/copper_fcc_lattice.hpp
new file mode 100644
index 0000000..e1b063b
--- /dev/null
+++ b/runtime/copper_fcc_lattice.hpp
@@ -0,0 +1,154 @@
+#include <iostream>
+#include <math.h>
+#include <mpi.h>
+//---
+#include "pairs.hpp"
+
+#pragma once
+
+namespace pairs {
+
+/* 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* seed) {
+ int k= (*seed) / IQ;
+ double ans;
+
+ *seed = IA * (*seed - k * IQ) - IR * k;
+ if(*seed < 0) *seed += IM;
+ ans = AM * (*seed);
+ return ans;
+}
+
+void random_reset(int *seed, int ibase, double *coord) {
+ int i;
+ char *str = (char *) &ibase;
+ int n = sizeof(int);
+ unsigned int hash = 0;
+
+ for (i = 0; i < n; i++) {
+ hash += str[i];
+ hash += (hash << 10);
+ hash ^= (hash >> 6);
+ }
+
+ str = (char *) coord;
+ n = 3 * sizeof(double);
+ for (i = 0; i < n; i++) {
+ hash += str[i];
+ hash += (hash << 10);
+ hash ^= (hash >> 6);
+ }
+
+ hash += (hash << 3);
+ hash ^= (hash >> 11);
+ hash += (hash << 15);
+
+ // keep 31 bits of unsigned int as new seed
+ // do not allow seed = 0, since will cause hang in gaussian()
+
+ *seed = hash & 0x7ffffff;
+ if (!(*seed)) *seed = 1;
+
+ // warm up the RNG
+
+ for (i = 0; i < 5; i++) myrandom(seed);
+ //save = 0;
+}
+
+double copper_fcc_lattice(PairsSimulation *ps, int nx, int ny, int nz, double xprd, double yprd, double zprd, double rho, int ntypes) {
+ auto shape = ps->getAsIntegerProperty(ps->getPropertyByName("shape"));
+ auto types = ps->getAsIntegerProperty(ps->getPropertyByName("type"));
+ auto flags = ps->getAsIntegerProperty(ps->getPropertyByName("flags"));
+ auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
+ auto positions = ps->getAsVectorProperty(ps->getPropertyByName("position"));
+ auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
+ double xlo = 0.0, xhi = xprd;
+ double ylo = 0.0, yhi = yprd;
+ double zlo = 0.0, zhi = zprd;
+ int natoms = 0;
+ int natoms_expected = 4 * nx * ny * nz;
+
+ double alat = pow((4.0 / 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 * nx - 1);
+ jlo = MAX(jlo, 0);
+ jhi = MIN(jhi, 2 * ny - 1);
+ klo = MAX(klo, 0);
+ khi = MIN(khi, 2 * nz - 1);
+
+ double 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(ps->getDomainPartitioner()->isWithinSubdomain(xtmp, ytmp, ztmp)) {
+ n = k * (2 * ny) * (2 * nx) + j * (2 * 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);
+
+ masses(natoms) = 1.0;
+ positions(natoms, 0) = xtmp;
+ positions(natoms, 1) = ytmp;
+ positions(natoms, 2) = ztmp;
+ velocities(natoms, 0) = vxtmp;
+ velocities(natoms, 1) = vytmp;
+ velocities(natoms, 2) = vztmp;
+ types(natoms) = rand() % ntypes;
+ flags(natoms) = 0;
+ shape(natoms) = 2; // point mass
+ natoms++;
+ }
+ }
+
+ 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++; }
+ }
+
+ if(natoms != natoms_expected) {
+ std::cerr << "copper_fcc_lattice(): incorrect number of atoms "
+ << "(" << natoms << " / " << natoms_expected << ")" << std::endl;
+ exit(-1);
+ }
+
+ return natoms;
+}
+
+}
diff --git a/runtime/thermo.hpp b/runtime/thermo.hpp
index 4afec58..7c999d6 100644
--- a/runtime/thermo.hpp
+++ b/runtime/thermo.hpp
@@ -1,4 +1,5 @@
#include <iostream>
+#include <math.h>
#include <mpi.h>
//---
#include "pairs.hpp"
@@ -7,13 +8,17 @@
namespace pairs {
-void compute_thermo(PairsSimulation *ps, int nlocal) {
+double compute_thermo(PairsSimulation *ps, int nlocal, double xprd, double yprd, double zprd, int print) {
auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
- const int natoms = 131072;
- const double xprd = 53.747078;
- const double yprd = 53.747078;
- const double zprd = 53.747078;
+ int natoms = nlocal;
+
+ if(ps->getDomainPartitioner()->getWorldSize() > 1) {
+ int global_natoms;
+ MPI_Reduce(&natoms, &global_natoms, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
+ natoms = global_natoms;
+ }
+
const double mvv2e = 1.0;
const double dof_boltz = (natoms * 3 - 3);
const double t_scale = mvv2e / dof_boltz;
@@ -36,11 +41,59 @@ void compute_thermo(PairsSimulation *ps, int nlocal) {
t = global_t;
}
- if(ps->getDomainPartitioner()->getRank() == 0) {
- t = t * t_scale;
+ t = t * t_scale;
+ if(print == 1 && ps->getDomainPartitioner()->getRank() == 0) {
p = (t * dof_boltz) * p_scale;
std::cout << t << "\t" << p << std::endl;
}
+
+ return t;
+}
+
+void adjust_thermo(PairsSimulation *ps, int nlocal, double xprd, double yprd, double zprd, double temp) {
+ auto velocities = ps->getAsVectorProperty(ps->getPropertyByName("linear_velocity"));
+ double vxtot = 0.0;
+ double vytot = 0.0;
+ double vztot = 0.0;
+ double tmp;
+ int natoms = nlocal;
+
+ for(int i = 0; i < nlocal; i++) {
+ vxtot += velocities(i, 0);
+ vytot += velocities(i, 1);
+ vztot += velocities(i, 2);
+ }
+
+ if(ps->getDomainPartitioner()->getWorldSize() > 1) {
+ int global_natoms;
+ MPI_Reduce(&natoms, &global_natoms, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
+ natoms = global_natoms;
+ MPI_Allreduce(&vxtot, &tmp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+ vxtot = tmp / natoms;
+ MPI_Allreduce(&vytot, &tmp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+ vytot = tmp / natoms;
+ MPI_Allreduce(&vztot, &tmp, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+ vztot = tmp / natoms;
+ } else {
+ vxtot /= natoms;
+ vytot /= natoms;
+ vztot /= natoms;
+ }
+
+ for(int i = 0; i < nlocal; i++) {
+ velocities(i, 0) -= vxtot;
+ velocities(i, 1) -= vytot;
+ velocities(i, 2) -= vztot;
+ }
+
+ double t = pairs::compute_thermo(ps, nlocal, xprd, yprd, zprd, 0);
+ double factor = sqrt(temp / t);
+
+ for(int i = 0; i < nlocal; i++) {
+ velocities(i, 0) *= factor;
+ velocities(i, 1) *= factor;
+ velocities(i, 2) *= factor;
+ }
}
}
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index 5a17c30..dba78f3 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -67,6 +67,7 @@ class CGen:
self.print("#include <stdlib.h>")
self.print("//---")
self.print("#include \"runtime/likwid-marker.h\"")
+ self.print("#include \"runtime/copper_fcc_lattice.hpp\"")
self.print("#include \"runtime/pairs.hpp\"")
self.print("#include \"runtime/read_from_file.hpp\"")
self.print("#include \"runtime/timing.hpp\"")
diff --git a/src/pairs/sim/copper_fcc_lattice.py b/src/pairs/sim/copper_fcc_lattice.py
new file mode 100644
index 0000000..ba966a5
--- /dev/null
+++ b/src/pairs/sim/copper_fcc_lattice.py
@@ -0,0 +1,34 @@
+from pairs.ir.assign import Assign
+from pairs.ir.block import pairs_inline
+from pairs.ir.functions import Call_Int, Call_Void
+from pairs.ir.particle_attributes import ParticleAttributeList
+from pairs.ir.types import Types
+from pairs.sim.grid import Grid3D
+from pairs.sim.lowerable import Lowerable
+
+
+class CopperFCCLattice(Lowerable):
+ def __init__(self, sim, nx, ny, nz, rho, temperature, ntypes):
+ super().__init__(sim)
+ self._nx = nx
+ self._ny = ny
+ self._nz = nz
+ self._rho = rho
+ self._temperature = temperature
+ self._ntypes = ntypes
+ lattice = pow((4.0 / rho), (1.0 / 3.0))
+ self._xprd = nx * lattice
+ self._yprd = ny * lattice
+ self._zprd = nz * lattice
+ sim.set_domain([0.0, 0.0, 0.0, self._xprd, self._yprd, self._zprd])
+
+ @pairs_inline
+ def lower(self):
+ Assign(self.sim, self.sim.nlocal,
+ Call_Int(self.sim, "pairs::copper_fcc_lattice",
+ [self._nx, self._ny, self._nz,
+ self._xprd, self._yprd, self._zprd,
+ self._rho, self._ntypes]))
+
+ Call_Void(self.sim, "pairs::adjust_thermo",
+ [self.sim.nlocal, self._xprd, self._yprd, self._zprd, self._temperature])
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 2ef4fea..c32c935 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -15,6 +15,7 @@ from pairs.sim.arrays import DeclareArrays
from pairs.sim.cell_lists import CellLists, BuildCellLists, BuildCellListsStencil, PartitionCellLists
from pairs.sim.comm import Comm
from pairs.sim.contact_history import ContactHistory, BuildContactHistory, ClearUnusedContactHistory, ResetContactHistoryUsageStatus
+from pairs.sim.copper_fcc_lattice import CopperFCCLattice
from pairs.sim.domain import InitializeDomain
from pairs.sim.domain_partitioning import DimensionRanges
from pairs.sim.features import AllocateFeatureProperties
@@ -24,6 +25,7 @@ from pairs.sim.lattice import ParticleLattice
from pairs.sim.neighbor_lists import NeighborLists, BuildNeighborLists
from pairs.sim.properties import AllocateProperties, AllocateContactProperties, ResetVolatileProperties
from pairs.sim.read_from_file import ReadParticleData
+from pairs.sim.thermo import ComputeThermo
from pairs.sim.timestep import Timestep
from pairs.sim.variables import DeclareVariables
from pairs.sim.vtk import VTKWrite
@@ -83,6 +85,7 @@ class Simulation:
self._compute_half = False
self._apply_list = None
self._enable_profiler = False
+ self._compute_thermo = 0
def enable_profiler(self):
self._enable_profiler = True
@@ -210,6 +213,9 @@ class Simulation:
props = [self.property(prop_name) for prop_name in prop_names]
self.setups.add_statement(ReadParticleData(self, filename, props, shape_id))
+ def copper_fcc_lattice(self, nx, ny, nz, rho, temperature, ntypes):
+ self.setups.add_statement(CopperFCCLattice(self, nx, ny, nz, rho, temperature, ntypes))
+
def build_cell_lists(self, spacing):
self.cell_lists = CellLists(self, self.grid, spacing, spacing)
return self.cell_lists
@@ -335,6 +341,9 @@ class Simulation:
def domain_partitioning(self):
return self._dom_part
+ def compute_thermo(self, every=0):
+ self._compute_thermo = every
+
def generate(self):
assert self._target is not None, "Target not specified!"
comm = Comm(self, self._dom_part)
@@ -354,7 +363,8 @@ class Simulation:
if self._use_contact_history:
if self.neighbor_lists is not None:
timestep_procedures.append(
- (BuildContactHistory(self, self._contact_history, self.cell_lists), every_reneighbor_params))
+ (BuildContactHistory(self, self._contact_history, self.cell_lists),
+ every_reneighbor_params))
timestep_procedures.append(ResetContactHistoryUsageStatus(self, self._contact_history))
@@ -363,6 +373,10 @@ class Simulation:
if self._use_contact_history:
timestep_procedures.append(ClearUnusedContactHistory(self, self._contact_history))
+ if self._compute_thermo != 0:
+ timestep_procedures.append(
+ (ComputeThermo(self), {'every': self._compute_thermo}))
+
timestep = Timestep(self, self.ntimesteps, timestep_procedures)
self.enter(timestep.block)
diff --git a/src/pairs/sim/thermo.py b/src/pairs/sim/thermo.py
new file mode 100644
index 0000000..557a756
--- /dev/null
+++ b/src/pairs/sim/thermo.py
@@ -0,0 +1,19 @@
+from pairs.ir.assign import Assign
+from pairs.ir.block import pairs_inline
+from pairs.ir.functions import Call_Int, Call_Void
+from pairs.ir.particle_attributes import ParticleAttributeList
+from pairs.ir.types import Types
+from pairs.sim.grid import Grid3D
+from pairs.sim.lowerable import Lowerable
+
+
+class ComputeThermo(Lowerable):
+ def __init__(self, sim):
+ super().__init__(sim)
+
+ @pairs_inline
+ def lower(self):
+ xprd = self.sim.grid.length(0)
+ yprd = self.sim.grid.length(1)
+ zprd = self.sim.grid.length(2)
+ Call_Void(self.sim, "pairs::compute_thermo", [self.sim.nlocal, xprd, yprd, zprd, 1])
--
GitLab