From 1eb665b4d61880eeb8eedd2cc4fbd91f4ceeaae1 Mon Sep 17 00:00:00 2001
From: Rafael Ravedutti <rafaelravedutti@gmail.com>
Date: Thu, 5 Oct 2023 23:13:21 +0200
Subject: [PATCH] Update data files and initialize domain separately from
 particle data

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
---
 data/minimd_setup_32x32x32.input         |  1 -
 data/minimd_setup_32x32x32_onetype.input |  1 -
 data/minimd_setup_4x4x4_onetype.input    |  1 -
 examples/dem.py                          | 75 ++++++++++++++++++++----
 examples/lj.py                           |  2 +
 runtime/read_from_file.hpp               |  2 +-
 src/pairs/sim/domain.py                  | 17 ++++++
 src/pairs/sim/grid.py                    | 19 +++---
 src/pairs/sim/read_from_file.py          | 11 ----
 src/pairs/sim/simulation.py              | 19 ++----
 10 files changed, 102 insertions(+), 46 deletions(-)
 create mode 100644 src/pairs/sim/domain.py

diff --git a/data/minimd_setup_32x32x32.input b/data/minimd_setup_32x32x32.input
index 72b1bc3..f8caf12 100644
--- a/data/minimd_setup_32x32x32.input
+++ b/data/minimd_setup_32x32x32.input
@@ -1,4 +1,3 @@
-0,53.747078,0,53.747078,0,53.747078
 2,1,0.000000,0.000000,0.000000,-1.168438,0.080847,-2.046589,0
 2,1,1.679596,0.000000,0.000000,0.651856,0.242263,-1.982591,0
 0,1,3.359192,0.000000,0.000000,-1.684546,0.403678,-1.918593,0
diff --git a/data/minimd_setup_32x32x32_onetype.input b/data/minimd_setup_32x32x32_onetype.input
index fb0d7be..80550bd 100644
--- a/data/minimd_setup_32x32x32_onetype.input
+++ b/data/minimd_setup_32x32x32_onetype.input
@@ -1,4 +1,3 @@
-0,53.747078,0,53.747078,0,53.747078
 1,0.000000,0.000000,0.000000,-1.168438,0.080847,-2.046589
 1,1.679596,0.000000,0.000000,0.651856,0.242263,-1.982591
 1,3.359192,0.000000,0.000000,-1.684546,0.403678,-1.918593
diff --git a/data/minimd_setup_4x4x4_onetype.input b/data/minimd_setup_4x4x4_onetype.input
index 239cf84..65c1672 100644
--- a/data/minimd_setup_4x4x4_onetype.input
+++ b/data/minimd_setup_4x4x4_onetype.input
@@ -1,4 +1,3 @@
-0,6.7184,0,6.7184,0,6.7184
 1,0,0,0,-1.13216,0.0489689,-2.03943
 1,1.6796,0,0,0.689294,0.210488,-1.97539
 1,3.3592,0,0,-1.6486,0.372007,-1.91135
diff --git a/examples/dem.py b/examples/dem.py
index 7d08a0f..e3d6a30 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -68,9 +68,51 @@ if target != 'cpu' and target != 'gpu':
     print(f"Invalid target, use {cmd} <cpu/gpu>")
 
 
-dt = 0.005
-cutoff_radius = 2.5
-skin = 0.3
+# BedGeneration {
+#    domainSize_SI < 0.8, 0.015, 0.2 >;
+#    blocks < 3, 3, 1 >;
+#    diameter_SI 0.0029;
+#    gravity_SI 9.81;
+#    densityFluid_SI 1000;
+#    densityParticle_SI 2550;
+#    generationSpacing_SI 0.005;
+#    initialVelocity_SI 1;
+#    dt_SI 5e-5;
+#    frictionCoefficient 0.5;
+#    restitutionCoefficient 0.1;
+#    collisionTime_SI 5e-4;
+#    poissonsRatio 0.22;
+#    timeSteps 10000;
+#    visSpacing 100;
+#    outFileName spheres_out.dat;
+#    denseBottomLayer False;
+#    bottomLayerOffsetFactor 1.0;
+#}
+
+# Config file parameters
+domainSize_SI = [0.8, 0.015, 0.2]
+blocks = [3, 3, 1]
+diameter_SI = 0.0029
+gravity_SI = 9.81
+densityFluid_SI = 1000
+densityParticle_SI = 2550
+generationSpacing_SI = 0.005
+initialVelocity_SI = 1
+dt_SI = 5e-5
+frictionCoefficient = 0.5
+restitutionCoefficient = 0.1
+collisionTime_SI = 5e-4
+poissonsRatio = 0.22
+timeSteps = 10000
+visSpacing = 100
+denseBottomLayer = False
+bottomLayerOffsetFactor = 1.0
+
+minDiameter_SI = diameter_SI * 0.9
+maxDiameter_SI = diameter_SI * 1.1
+linkedCellWidth = 1.01 * maxDiameter_SI
+
+skin = 0.1
 ntypes = 1
 stiffness_norm = 1.0
 stiffness_tan = 1.0
@@ -79,10 +121,6 @@ damping_tan = 1.0
 friction_static = 1.0
 friction_dynamic = 1.0
 
-densityParticle_SI = 1.0
-densityFluid_SI = 1.0
-gravity_SI = 1.0
-
 psim = pairs.simulation("dem", debug=True)
 psim.add_position('position')
 psim.add_property('mass', pairs.double(), 1.0)
@@ -102,11 +140,26 @@ psim.add_contact_property('is_sticking', pairs.int32(), 0)
 psim.add_contact_property('tangential_spring_displacement', pairs.vector(), [0.0, 0.0, 0.0])
 psim.add_contact_property('impact_velocity_magnitude', pairs.double(), 0.0)
 
-psim.read_particle_data("data/fluidized_bed.input", ['mass', 'position', 'linear_velocity'], pairs.sphere())
-psim.build_neighbor_lists(cutoff_radius + skin)
+psim.set_domain([0.0, 0.0, 0.0, domainSize_SI[0], domainSize_SI[1], domainSize_SI[2]])
+psim.read_particle_data(
+    "data/spheres.input",
+    ['type', 'mass', 'radius', 'position', 'linear_velocity', 'flags'],
+    pairs.sphere())
+
+psim.read_particle_data(
+    "data/spheres_bottom.input",
+    ['type', 'mass', 'radius', 'position', 'linear_velocity', 'flags'],
+    pairs.sphere())
+
+psim.read_particle_data(
+    "data/planes.input",
+    ['type', 'mass', 'position', 'normal', 'flags'],
+    pairs.halfspace())
+
+psim.build_neighbor_lists(linkedCellWidth + skin)
 psim.vtk_output(f"output/test_{target}")
-psim.compute(linear_spring_dashpot, cutoff_radius, symbols={'dt': dt})
-psim.compute(euler, symbols={'dt': dt})
+psim.compute(linear_spring_dashpot, linkedCellWidth + skin, symbols={'dt': dt_SI})
+psim.compute(euler, symbols={'dt': dt_SI})
 psim.compute(gravity, symbols={'densityParticle_SI': densityParticle_SI,
                                'densityFluid_SI': densityFluid_SI,
                                'gravity_SI': gravity_SI,
diff --git a/examples/lj.py b/examples/lj.py
index 4fabc8a..b280027 100644
--- a/examples/lj.py
+++ b/examples/lj.py
@@ -34,6 +34,8 @@ psim.add_property('force', pairs.vector(), volatile=True)
 psim.add_feature('type', ntypes)
 psim.add_feature_property('type', 'epsilon', pairs.double(), [sigma for i in range(ntypes * ntypes)])
 psim.add_feature_property('type', 'sigma6', pairs.double(), [epsilon for i in range(ntypes * ntypes)])
+#psim.set_domain(0.0, 0.0, 0.0, 6.7184, 6.7184, 6.7184) # for 4x4x4 setup
+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.sphere())
 psim.build_neighbor_lists(cutoff_radius + skin)
 psim.vtk_output(f"output/test_{target}")
diff --git a/runtime/read_from_file.hpp b/runtime/read_from_file.hpp
index 95f007e..e29e674 100644
--- a/runtime/read_from_file.hpp
+++ b/runtime/read_from_file.hpp
@@ -37,7 +37,7 @@ size_t read_particle_data(PairsSimulation *ps, const char *filename, const prope
     size_t n = 0;
 
     if(in_file.is_open()) {
-        std::getline(in_file, line);
+        //std::getline(in_file, line);
         while(std::getline(in_file, line)) {
             std::stringstream line_stream(line);
             std::string in0;
diff --git a/src/pairs/sim/domain.py b/src/pairs/sim/domain.py
new file mode 100644
index 0000000..a306353
--- /dev/null
+++ b/src/pairs/sim/domain.py
@@ -0,0 +1,17 @@
+from pairs.ir.block import pairs_inline
+from pairs.ir.functions import Call_Void
+from pairs.ir.types import Types
+from pairs.sim.lowerable import Lowerable
+
+
+class InitializeDomain(Lowerable):
+    def __init__(self, sim):
+        super().__init__(sim)
+
+    @pairs_inline
+    def lower(self):
+        dom_part = self.sim.domain_partitioning()
+        grid_array = [(self.sim.grid.min(d), self.sim.grid.max(d)) for d in range(self.sim.ndims())]
+        Call_Void(self.sim, "pairs->initDomain", [param for delim in grid_array for param in delim]),
+        Call_Void(self.sim, "pairs->fillCommunicationArrays", [dom_part.neighbor_ranks, dom_part.pbc, dom_part.subdom])
+
diff --git a/src/pairs/sim/grid.py b/src/pairs/sim/grid.py
index 54c8e41..805ade2 100644
--- a/src/pairs/sim/grid.py
+++ b/src/pairs/sim/grid.py
@@ -1,3 +1,4 @@
+from pairs.ir.lit import Lit
 from pairs.ir.types import Types
 
 
@@ -5,7 +6,7 @@ class Grid:
     def __init__(self, sim, config):
         self.sim = sim
         self.ndims = len(config)
-        self.config = config
+        self.config = [(Lit.cvt(sim, dmin), Lit.cvt(sim, dmax)) for dmin, dmax in config]
 
     def range(self, dim):
         return self.config[dim]
@@ -39,14 +40,14 @@ class Grid:
 
 
 class Grid2D(Grid):
-    def __init__(self, sim, xmin, xmax, ymin, ymax):
-        config = [[xmin, xmax], [ymin, ymax]]
+    def __init__(self, sim, xmin, ymin, xmax, ymax):
+        config = [(xmin, xmax), (ymin, ymax)]
         super().__init__(sim, config)
 
 
 class Grid3D(Grid):
-    def __init__(self, sim, xmin, xmax, ymin, ymax, zmin, zmax):
-        config = [[xmin, xmax], [ymin, ymax], [zmin, zmax]]
+    def __init__(self, sim, xmin, ymin, zmin, xmax, ymax, zmax):
+        config = [(xmin, xmax), (ymin, ymax), (zmin, zmax)]
         super().__init__(sim, config)
 
 
@@ -56,6 +57,10 @@ class MutableGrid(Grid):
     def __init__(self, sim, ndims):
         self.id = MutableGrid.last_id
         prefix = f"grid{self.id}_"
-        config = [[sim.add_var(f"{prefix}d{d}_min", Types.Double), sim.add_var(f"{prefix}d{d}_max", Types.Double)] for d in range(ndims)]
-        super().__init__(sim, config)
+        super().__init__(sim, [
+            (sim.add_var(f"{prefix}d{d}_min", Types.Double),
+             sim.add_var(f"{prefix}d{d}_max", Types.Double))
+            for d in range(ndims)
+        ])
+
         MutableGrid.last_id += 1
diff --git a/src/pairs/sim/read_from_file.py b/src/pairs/sim/read_from_file.py
index c87171a..4c3815c 100644
--- a/src/pairs/sim/read_from_file.py
+++ b/src/pairs/sim/read_from_file.py
@@ -12,21 +12,10 @@ class ReadParticleData(Lowerable):
         super().__init__(sim)
         self.filename = filename
         self.attrs = ParticleAttributeList(sim, items)
-        self.grid = MutableGrid(sim, sim.ndims())
-        self.grid_buffer = sim.add_static_array("grid_buffer", [sim.ndims() * 2], Types.Double)
         self.shape_id = shape_id
 
     @pairs_inline
     def lower(self):
-        Call_Void(self.sim, "pairs::read_grid_data", [self.filename, self.grid_buffer])
-        for d in range(self.sim.ndims()):
-            Assign(self.sim, self.grid.min(d), self.grid_buffer[d * 2 + 0])
-            Assign(self.sim, self.grid.max(d), self.grid_buffer[d * 2 + 1])
-
-        dom_part = self.sim.domain_partitioning()
-        grid_array = [[self.grid.min(d), self.grid.max(d)] for d in range(self.sim.ndims())]
-        Call_Void(self.sim, "pairs->initDomain", [param for delim in grid_array for param in delim]),
-        Call_Void(self.sim, "pairs->fillCommunicationArrays", [dom_part.neighbor_ranks, dom_part.pbc, dom_part.subdom])
         Assign(self.sim, self.sim.nlocal, Call_Int(self.sim, "pairs::read_particle_data", [self.filename, self.attrs, self.attrs.length(), self.shape_id]))
 
 
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 328b961..edc321f 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -16,6 +16,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
+from pairs.sim.domain import InitializeDomain
 from pairs.sim.domain_partitioning import DimensionRanges
 from pairs.sim.features import AllocateFeatureProperties
 from pairs.sim.grid import Grid2D, Grid3D
@@ -167,24 +168,16 @@ class Simulation:
     def var(self, var_name):
         return self.vars.find(var_name)
 
-    def grid_2d(self, xmin, xmax, ymin, ymax):
-        self.grid = Grid2D(self, xmin, xmax, ymin, ymax)
-        return self.grid
-
-    def grid_3d(self, xmin, xmax, ymin, ymax, zmin, zmax):
-        self.grid = Grid3D(self, xmin, xmax, ymin, ymax, zmin, zmax)
-        return self.grid
+    def set_domain(self, grid):
+        self.grid = Grid3D(self, grid[0], grid[1], grid[2], grid[3], grid[4], grid[5])
+        self.setups.add_statement(InitializeDomain(self))
 
     def create_particle_lattice(self, grid, spacing, props={}):
-        positions = self.position()
-        lattice = ParticleLattice(self, grid, spacing, props, positions)
-        self.setups.add_statement(lattice)
+        self.setups.add_statement(ParticleLattice(self, grid, spacing, props, self.position()))
 
     def read_particle_data(self, filename, prop_names, shape_id):
         props = [self.property(prop_name) for prop_name in prop_names]
-        read_object = ReadParticleData(self, filename, props, shape_id)
-        self.setups.add_statement(read_object)
-        self.grid = read_object.grid
+        self.setups.add_statement(ReadParticleData(self, filename, props, shape_id))
 
     def build_cell_lists(self, spacing):
         self.cell_lists = CellLists(self, self.grid, spacing, spacing)
-- 
GitLab