diff --git a/dem/app/app.py b/dem/app/app.py index a4a6f73..68c15b4 100644 --- a/dem/app/app.py +++ b/dem/app/app.py @@ -8,6 +8,7 @@ from dem.app.silo import * from dem.app.mill import * from dem.app.circular import * +from dem.app.silo_open import * # Declare factory app_factory = factory() @@ -21,3 +22,4 @@ app_factory.register("silo", silo) app_factory.register("mill", mill) app_factory.register("circular", circular) +app_factory.register("silo_open", silo_open) diff --git a/dem/app/base_app.py b/dem/app/base_app.py index 306fad6..64919c3 100644 --- a/dem/app/base_app.py +++ b/dem/app/base_app.py @@ -51,6 +51,18 @@ def update(self): self.it += 1 self.t += self.dt + ### ************************************************ + ### Check if particles are still in the domain + def check_particles(self, d): + + lst = [] + + for i in range(self.p.np): + if (not d.is_in(self.p.x[i,:])): + lst.append(i) + + self.p.delete(lst) + ### ************************************************ ### Plot def plot(self): diff --git a/dem/app/silo_open.py b/dem/app/silo_open.py new file mode 100644 index 0000000..a2527a3 --- /dev/null +++ b/dem/app/silo_open.py @@ -0,0 +1,127 @@ +# Generic imports +import os +import math +import random + +# Custom imports +from dem.app.base_app import * + +### ************************************************ +### Silo case with open boundary +class silo_open(base_app): + ### ************************************************ + ### Constructor + def __init__(self, + name = 'silo_open', + t_max = 3.5, + dt = 2.5e-5, + angle = 0.0, + plot_freq = 500, + plot_show = False, + plot_trajectory = False, + plot_png = True): + super().__init__() + + self.name = name + self.t_max = t_max + self.dt = dt + self.plot_freq = plot_freq + self.plot_show = plot_show + self.plot_trajectory = plot_trajectory + self.plot_png = plot_png + + self.nt = int(self.t_max/self.dt) + self.plot_it = 0 + self.check_freq = 10 # check particles every 10 iterations + + self.n_row = 25 # nb of particles on a row at start + self.n_col = 20 # nb of particles on a col at start + self.radius = 0.025 + + self.p = particles(np = self.n_row*self.n_col, + nt = self.nt, + material = "steel", + radius = self.radius, + color = "b", + store = False, + search = "nearest", + rad_coeff = 2.0) + + self.p.e_wall[:] = 0.5 + self.p.e_part[:] = 0.5 + + colors = np.array(['r', 'g', 'b', 'c', 'm', 'y', 'k']) + self.p.c = colors[np.random.randint(0,len(colors),size=self.p.np)] + + self.d = domain_factory.create("rectangle", + x_min = 0.0, + x_max = 3.0, + y_min = 0.0, + y_max = 5.0, + angle = 0.0, + material = "steel", + open_boundary = [True, False, False, False]) + + self.o0 = domain_factory.create("rectangle", + x_min =-1.0, + x_max = 1.4, + y_min = 2.5, + y_max = 2.6, + angle =-30.0, + plot_fill = True, + material = "steel") + self.o1 = domain_factory.create("rectangle", + x_min = 1.6, + x_max = 4.0, + y_min = 2.5, + y_max = 2.6, + angle = 30.0, + plot_fill = True, + material = "steel") + + self.d_lst = [self.d, self.o0, self.o1] + + self.path = self.base_path+'/'+self.name + os.makedirs(self.path, exist_ok=True) + + self.reset() + + ### ************************************************ + ### Reset app + def reset(self): + + self.it = 0 + self.t = 0.0 + sep = 4.0*self.radius + mid = 0.5*(self.d.x_min + self.d.x_max) + half = 0.5*self.n_row*(self.radius + 0.75*sep) + + for i in range(self.n_row): + for j in range(self.n_col): + self.p.x[self.n_col*i+j,0] = mid - half + sep*i + self.radius*random.random() + self.p.x[self.n_col*i+j,1] = self.d.y_max - 2*sep - sep*j + + ### ************************************************ + ### Compute forces + def forces(self): + + # Removing lost particles requires to recompute the + # nearest neighbor lists. Everytime the check is done, + # we force the recomputation + force_nearest = False + if (self.it%self.check_freq == 0): + self.check_particles(self.d) + force_nearest = True + + self.p.reset_forces() + self.p.collisions(self.dt, force_nearest=force_nearest) + self.d.collisions(self.p, self.dt) + self.o0.collisions(self.p, self.dt) + self.o1.collisions(self.p, self.dt) + self.p.gravity(self.g) + + ### ************************************************ + ### Iteration printings (overrides base_app::printings) + def printings(self): + + print("# it = "+str(self.it)+", t = {0:.3f}".format(self.t)+" / {0:.3f}".format(self.t_max)+", n_particles = "+str(self.p.np), end='\r') diff --git a/dem/src/core/particles.py b/dem/src/core/particles.py index 627965e..a738c7d 100644 --- a/dem/src/core/particles.py +++ b/dem/src/core/particles.py @@ -86,17 +86,43 @@ def set_material(self, i, m): ### Compute maximal velocity def max_velocity(self): - return np.max(np.linalg.norm(self.v, axis=1)) + if (len(self.v) == 0): + return 0.0 + else: + return np.max(np.linalg.norm(self.v, axis=1)) ### ************************************************ ### Compute maximal radius def max_radius(self): - return np.max(self.r) + if (len(self.r) == 0): + return 0.0 + else: + return np.max(self.r) + + ### ************************************************ + ### Remove a list of particles + def delete(self, lst): + + self.np -= len(lst) + self.m = np.delete(self.m, lst, axis=0) + self.r = np.delete(self.r, lst, axis=0) + self.x = np.delete(self.x, lst, axis=0) + self.d = np.delete(self.d, lst, axis=0) + self.v = np.delete(self.v, lst, axis=0) + self.a = np.delete(self.a, lst, axis=0) + self.f = np.delete(self.f, lst, axis=0) + self.e_wall = np.delete(self.e_wall, lst, axis=0) + self.mu_wall = np.delete(self.mu_wall, lst, axis=0) + self.e_part = np.delete(self.e_part, lst, axis=0) + self.mu_part = np.delete(self.mu_part, lst, axis=0) + self.Y = np.delete(self.Y, lst, axis=0) + self.G = np.delete(self.G, lst, axis=0) + self.c = np.delete(np.array(self.c), lst, axis=0).tolist() ### ************************************************ ### Compute collisions between particles - def collisions(self, dt): + def collisions(self, dt, force_nearest=False): if (self.search == "linear"): @@ -117,7 +143,7 @@ def collisions(self, dt): # If list is not valid anymore v_max = self.max_velocity() self.d_ngb += 2.0*v_max*dt - if (self.d_ngb > 0.99*self.m_rad): + if ((self.d_ngb > 0.99*self.m_rad) or force_nearest): self.d_ngb = 0.0 self.ngbi, self.ngbj = list_ngbs(self.np, self.x, self.r, self.m_rad) diff --git a/dem/src/domain/base.py b/dem/src/domain/base.py index c9457b8..56cbf59 100644 --- a/dem/src/domain/base.py +++ b/dem/src/domain/base.py @@ -9,19 +9,7 @@ def __init__(self): pass ### ************************************************ - ### Compute distance to particle - def distance(self, p): - - raise NotImplementedError - - ### ************************************************ - ### Detect collision with particle - def collide(self, p): - - raise NotImplementedError - - ### ************************************************ - ### Plot domain - def plot(self, ax): + ### Check if point is in domain + def is_in(self, pt): raise NotImplementedError diff --git a/dem/src/domain/rectangle.py b/dem/src/domain/rectangle.py index 579f0dc..46abd02 100644 --- a/dem/src/domain/rectangle.py +++ b/dem/src/domain/rectangle.py @@ -17,13 +17,14 @@ class rectangle(base_domain): ### ************************************************ ### Constructor def __init__(self, - x_min = 0.0, - x_max = 1.0, - y_min = 0.0, - y_max = 1.0, - angle = 0.0, # in degrees - plot_fill = False, - material = "steel"): + x_min = 0.0, + x_max = 1.0, + y_min = 0.0, + y_max = 1.0, + angle = 0.0, # in degrees + plot_fill = False, + material = "steel", + open_boundary = [False, False, False, False]): # Type self.type = "rectangle" @@ -44,6 +45,9 @@ def __init__(self, # Material self.mat = material_factory.create(material) + # Open boundaries + self.open_boundary = np.array(open_boundary) + # Ficticious parameters for domain self.r = 1.0e8 self.m = 1.0e8 @@ -81,7 +85,7 @@ def collisions(self, p, dt): # Search for collisions linearly and compute forces linear_search(self.seg_pts, self.r, self.m, - self.v, self.mat.Y, self.mat.G, + self.v, self.mat.Y, self.mat.G, self.open_boundary, p.x, p.r, p.m, p.v, p.e_wall, p.mu_wall, p.Y, p.G, p.f, p.np, dt) @@ -96,13 +100,34 @@ def rotate(self, p, pc, angle): dp[1] = (p[0]-pc[0])*sint + (p[1]-pc[1])*cost + pc[1] p[:] = dp[:] + ### ************************************************ + ### Check if point is in domain + ### pt is assumed to be an np array of size 2 + def is_in(self, pm): + + p1p2 = self.p2 - self.p1 + p1pm = pm - self.p1 + p1p4 = self.p4 - self.p1 + + p1p2p1pm = np.dot(p1p2, p1pm) + p1p2p1p2 = np.dot(p1p2, p1p2) + p1p4p1pm = np.dot(p1p4, p1pm) + p1p4p1p4 = np.dot(p1p4, p1p4) + + if ((p1p2p1pm > 0.0) and + (p1p2p1p2 > p1p2p1pm) and + (p1p4p1pm > 0.0) and + (p1p4p1p4 > p1p4p1pm)): return True + + return False + ### ************************************************ ### Distance from rectangle domain to given coordinates ### Prefix d_ corresponds to domain ### Prefix p_ corresponds to particle @nb.njit(cache=True) def linear_search(d_pts, d_r, d_m, - d_v, d_mat_Y, d_mat_G, + d_v, d_mat_Y, d_mat_G, open_boundary, p_x, p_r, p_m, p_v, p_e_wall, p_mu_wall, p_Y, p_G, p_f, n, dt): @@ -111,6 +136,10 @@ def linear_search(d_pts, d_r, d_m, # Loop on rectangle sides for j in range(4): + + # If this boundary is open, loop + if (open_boundary[j]): continue + dx, nrm = p_to_segment(p_x[i], d_pts[j,0], d_pts[j,1]) dx = dx - p_r[i]