diff --git a/README.rst b/README.rst index c6b2cc61d..206c308d7 100644 --- a/README.rst +++ b/README.rst @@ -94,6 +94,8 @@ basic examples. Some of the formulations available are: pp. 124--149) - `Akinci `_ (Akinci et al. ACM Trans. Graph., 2012, pp. 62:1--62:8) +- `Direct Fiber Simulation `_ for + fiber suspensions (Meyer et al. 2020, JCS, 4(2), pp. 77) Boundary conditions from the following papers are implemented: diff --git a/pysph/base/gpu_domain_manager.py b/pysph/base/gpu_domain_manager.py index e63449187..a395118d7 100644 --- a/pysph/base/gpu_domain_manager.py +++ b/pysph/base/gpu_domain_manager.py @@ -14,8 +14,10 @@ class GPUDomainManager(DomainManagerBase): def __init__(self, xmin=-1000., xmax=1000., ymin=0., ymax=0., zmin=0., zmax=0., - periodic_in_x=False, periodic_in_y=False, - periodic_in_z=False, n_layers=2.0, backend=None, props=None, + periodic_in_x=False, periodic_in_y=False, periodic_in_z=False, + gamma_yx=0.0, gamma_zx=0.0, gamma_zy=0.0, + n_layers=2.0, backend=None, props=None, + dt=0.0, calls_per_step=2, mirror_in_x=False, mirror_in_y=False, mirror_in_z=False): """Constructor""" DomainManagerBase.__init__(self, xmin=xmin, xmax=xmax, diff --git a/pysph/base/nnps_base.pxd b/pysph/base/nnps_base.pxd index 7e927194a..bb1af2aab 100644 --- a/pysph/base/nnps_base.pxd +++ b/pysph/base/nnps_base.pxd @@ -138,7 +138,7 @@ cpdef UIntArray arange_uint(int start, int stop=*) # Basic particle array wrapper used for NNPS cdef class NNPSParticleArrayWrapper: - cdef public DoubleArray x,y,z,h + cdef public DoubleArray x,y,z,u,v,w,h cdef public UIntArray gid cdef public IntArray tag cdef public ParticleArray pa @@ -164,6 +164,10 @@ cdef class DomainManagerBase: cdef public double ytranslate cdef public double ztranslate + cdef public double gamma_yx + cdef public double gamma_zx + cdef public double gamma_zy + cdef public int dim cdef public bint periodic_in_x, periodic_in_y, periodic_in_z cdef public bint is_periodic @@ -180,8 +184,14 @@ cdef class DomainManagerBase: cdef public double radius_scale # Radius scale for kernel cdef public double n_layers # Number of layers of ghost particles + cdef public double t # Time for Lees-Edwards BC + cdef public double dt # Time step for Lees-Edwards BC + cdef public int calls_per_step # Number of calls per step + cdef public int loops # loop counter + #cdef double dbl_max # Maximum value of double + # remove ghost particles from a previous iteration cpdef _remove_ghosts(self) @@ -199,6 +209,10 @@ cdef class CPUDomainManager(DomainManagerBase): # Convenience function to add a value to a carray cdef _add_to_array(self, DoubleArray arr, double disp, int start=*) + # Functions to shift ghost particle periodically in a direction + cdef _shift_periodic(self, DoubleArray arr, double disp, double min_pos, + double max_pos, int start=*) + # Convenience function to multiply a value to a carray cdef _mul_to_array(self, DoubleArray arr, double val) diff --git a/pysph/base/nnps_base.pyx b/pysph/base/nnps_base.pyx index 414d03a42..048b66cd0 100644 --- a/pysph/base/nnps_base.pyx +++ b/pysph/base/nnps_base.pyx @@ -225,6 +225,24 @@ cdef class NNPSParticleArrayWrapper: self.y = pa.get_carray('y') self.z = pa.get_carray('z') self.h = pa.get_carray('h') + if PyDict_Contains(pa.properties, 'u') == 1: + self.u = pa.get_carray('u') + else: + self.u = DoubleArray() + for i in range(self.x.length): + self.u.append(0.0) + if PyDict_Contains(pa.properties, 'v') == 1: + self.v = pa.get_carray('v') + else: + self.v = DoubleArray() + for i in range(self.y.length): + self.v.append(0.0) + if PyDict_Contains(pa.properties, 'w') == 1: + self.w = pa.get_carray('w') + else: + self.w = DoubleArray() + for i in range(self.z.length): + self.w.append(0.0) self.gid = pa.get_carray('gid') self.tag = pa.get_carray('tag') @@ -243,10 +261,11 @@ cdef class NNPSParticleArrayWrapper: cdef class DomainManager: def __init__(self, double xmin=-1000, double xmax=1000, double ymin=0, double ymax=0, double zmin=0, double zmax=0, + double gamma_yx=0.0, double gamma_zx=0.0,double gamma_zy=0.0, periodic_in_x=False, periodic_in_y=False, periodic_in_z=False, - double n_layers=2.0, backend=None, props=None, + double n_layers=2.0, backend=None, double dt=0.0, + int calls_per_step=2, props=None, mirror_in_x=False, mirror_in_y=False, mirror_in_z=False): - """Constructor Parameters @@ -279,9 +298,13 @@ cdef class DomainManager: domain_manager = CPUDomainManager self.manager = domain_manager( xmin=xmin, xmax=xmax, ymin=ymin, - ymax=ymax, zmin=zmin, zmax=zmax, periodic_in_x=periodic_in_x, - periodic_in_y=periodic_in_y, periodic_in_z=periodic_in_z, + ymax=ymax, zmin=zmin, zmax=zmax, + periodic_in_x=periodic_in_x, + periodic_in_y=periodic_in_y, + periodic_in_z=periodic_in_z, + gamma_yx=gamma_yx, gamma_zx=gamma_zx,gamma_zy=gamma_zy, n_layers=n_layers, backend=self.backend, props=props, + dt=dt, calls_per_step=calls_per_step, mirror_in_x=mirror_in_x, mirror_in_y=mirror_in_y, mirror_in_z=mirror_in_z ) @@ -318,8 +341,10 @@ cdef class DomainManagerBase(object): def __init__(self, double xmin=-1000, double xmax=1000, double ymin=0, double ymax=0, double zmin=0, double zmax=0, periodic_in_x=False, periodic_in_y=False, periodic_in_z=False, - double n_layers=2.0, props=None, mirror_in_x=False, - mirror_in_y=False, mirror_in_z=False): + double gamma_yx=0.0, double gamma_zx=0.0, double gamma_zy=0.0, + double n_layers=2.0, backend=None, props=None, + double dt=0.0, int calls_per_step=2, + mirror_in_x=False, mirror_in_y=False, mirror_in_z=False): """Constructor The n_layers argument specifies the number of ghost layers as multiples @@ -356,6 +381,11 @@ cdef class DomainManagerBase(object): self.ytranslate = ymax - ymin self.ztranslate = zmax - zmin + # Shear rate for Lees-Edwards condition + self.gamma_yx = gamma_yx + self.gamma_zx = gamma_zx + self.gamma_zy = gamma_zy + # empty list of particle array wrappers for now self.pa_wrappers = [] self.narrays = 0 @@ -367,11 +397,18 @@ cdef class DomainManagerBase(object): # default DomainManager in_parallel is set to False self.in_parallel = False + # time for Lees-Edwards BC + self.t = 0 + self.dt = dt + self.calls_per_step = calls_per_step + self.loops = 0 + def _check_limits(self, xmin, xmax, ymin, ymax, zmin, zmax): """Sanity check on the limits""" if ( (xmax < xmin) or (ymax < ymin) or (zmax < zmin) ): raise ValueError("Invalid domain limits!") + #### Public protocol ################################################ def set_pa_wrappers(self, wrappers): self.pa_wrappers = wrappers self.narrays = len(wrappers) @@ -435,7 +472,9 @@ cdef class CPUDomainManager(DomainManagerBase): def __init__(self, double xmin=-1000, double xmax=1000, double ymin=0, double ymax=0, double zmin=0, double zmax=0, periodic_in_x=False, periodic_in_y=False, periodic_in_z=False, + double gamma_yx=0.0, double gamma_zx=0.0, double gamma_zy=0.0, double n_layers=2.0, backend=None, props=None, + double dt=0.0, int calls_per_step=2, mirror_in_x=False, mirror_in_y=False, mirror_in_z=False): """Constructor @@ -452,7 +491,9 @@ cdef class CPUDomainManager(DomainManagerBase): self, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, periodic_in_x=periodic_in_x, periodic_in_y=periodic_in_y, - periodic_in_z=periodic_in_z, n_layers=n_layers, props=props, + periodic_in_z=periodic_in_z, gamma_yx=gamma_yx, gamma_zx=gamma_zx, + gamma_zy=gamma_zy, n_layers=n_layers, props=props, dt=dt, + calls_per_step=calls_per_step, mirror_in_x=mirror_in_x, mirror_in_y=mirror_in_y, mirror_in_z=mirror_in_z ) @@ -498,6 +539,14 @@ cdef class CPUDomainManager(DomainManagerBase): # Update GPU. self._update_gpu() + # update total time (Designed for Transport velocity step) + # print("Loops: %d" % self.loops) + # print("calls_per_step : %d" % self.calls_per_step) + if self.loops > 0 and (self.loops % self.calls_per_step) == 0: + self.t = self.t + self.dt + # print(self.t) + self.loops =self.loops + 1 + #### Private protocol ############################################### cdef _add_to_array(self, DoubleArray arr, double disp, int start=0): cdef int i @@ -519,6 +568,21 @@ cdef class CPUDomainManager(DomainManagerBase): for i in range(arr.length): arr.data[i] *= val + + cdef _shift_periodic(self, DoubleArray arr, double disp, double min_pos, + double max_pos, int start=0): + cdef int i + cdef double L = max_pos-min_pos + for i in range(arr.length - start): + if disp > 0.0: + arr.data[start + i] += disp % L + else: + arr.data[start + i] -= (-disp) % L + if arr.data[start + i] < min_pos: + arr.data[start + i] += L + elif arr.data[start + i] > max_pos: + arr.data[start + i] -= L + cdef _create_ghosts_mirror(self): """Identify boundary particles and create images. @@ -712,6 +776,9 @@ cdef class CPUDomainManager(DomainManagerBase): added.tag[:] = Ghost pa.append_parray(added) + @cython.wraparound (False) #turn off negative indexing + @cython.boundscheck(False) #turn off bounds-checking + @cython.cdivision(True) cdef _box_wrap_periodic(self): """Box-wrap particles for periodicity @@ -747,22 +814,63 @@ cdef class CPUDomainManager(DomainManagerBase): # iterate over each array and mark for translation for pa_wrapper in self.pa_wrappers: x = pa_wrapper.x; y = pa_wrapper.y; z = pa_wrapper.z + u = pa_wrapper.u; v = pa_wrapper.v; w = pa_wrapper.w np = x.length # iterate over particles and box-wrap for i in range(np): - if periodic_in_x: - if x.data[i] < xmin : x.data[i] = x.data[i] + xtranslate - if x.data[i] > xmax : x.data[i] = x.data[i] - xtranslate + if x.data[i] < xmin : + x.data[i] = x.data[i] + xtranslate + if periodic_in_y: + y.data[i] = ( + y.data[i] + (xtranslate*self.gamma_yx*self.t) % ytranslate + ) + v.data[i] = v.data[i] + xtranslate*self.gamma_yx + if periodic_in_z: + z.data[i] = ( + z.data[i] + (xtranslate*self.gamma_zx*self.t) % ztranslate + ) + w.data[i] = w.data[i] + xtranslate*self.gamma_zx + if x.data[i] > xmax : + x.data[i] = x.data[i] - xtranslate + if periodic_in_y: + y.data[i] = ( + y.data[i] - (xtranslate*self.gamma_yx*self.t) % ytranslate + ) + v.data[i] = v.data[i] - xtranslate*self.gamma_yx + if periodic_in_z: + z.data[i] = ( + z.data[i] - (xtranslate*self.gamma_zx*self.t) % ztranslate + ) + w.data[i] = w.data[i] - xtranslate*self.gamma_zx if periodic_in_y: - if y.data[i] < ymin : y.data[i] = y.data[i] + ytranslate - if y.data[i] > ymax : y.data[i] = y.data[i] - ytranslate + if y.data[i] < ymin : + y.data[i] = y.data[i] + ytranslate + if periodic_in_z: + z.data[i] = ( + z.data[i] + (ytranslate*self.gamma_zy*self.t) % ztranslate + ) + w.data[i] = w.data[i] + ytranslate*self.gamma_zy + if y.data[i] > ymax : + y.data[i] = y.data[i] - ytranslate + if periodic_in_z: + z.data[i] = ( + z.data[i] - (ytranslate*self.gamma_zy*self.t) % ztranslate + ) + w.data[i] = w.data[i] - ytranslate*self.gamma_zy if periodic_in_z: - if z.data[i] < zmin : z.data[i] = z.data[i] + ztranslate - if z.data[i] > zmax : z.data[i] = z.data[i] - ztranslate + if z.data[i] < zmin : + z.data[i] = z.data[i] + ztranslate + if z.data[i] > zmax : + z.data[i] = z.data[i] - ztranslate + + def _check_limits(self, xmin, xmax, ymin, ymax, zmin, zmax): + """Sanity check on the limits""" + if ( (xmax < xmin) or (ymax < ymin) or (zmax < zmin) ): + raise ValueError("Invalid domain limits!") cdef _create_ghosts_periodic(self): """Identify boundary particles and create images. @@ -803,6 +911,7 @@ cdef class CPUDomainManager(DomainManagerBase): cdef ParticleArray pa, ghost_pa cdef DoubleArray x, y, z, xt, yt, zt cdef double xi, yi, zi + cdef double shift, L cdef int array_index, i, np, start # temporary indices for particles to be replicated @@ -831,10 +940,12 @@ cdef class CPUDomainManager(DomainManagerBase): x = pa_wrapper.x; y = pa_wrapper.y; z = pa_wrapper.z # reset the length of the arrays - x_low.reset(); x_high.reset(); y_high.reset(); y_low.reset() + x_low.reset(); x_high.reset() + y_low.reset(); y_high.reset() z_low.reset(); z_high.reset() np = x.length + for i in range(np): xi = x.data[i]; yi = y.data[i]; zi = z.data[i] @@ -862,6 +973,20 @@ cdef class CPUDomainManager(DomainManagerBase): ) self._add_to_array(ghost_pa.get_carray('x'), xtranslate, start=start) + if periodic_in_y and self.gamma_yx: + shift = xtranslate*self.gamma_yx + L = self.ymin + U = self.ymax + self._shift_periodic(ghost_pa.get_carray('y'), + shift*self.t, L, U, start=start) + self._add_to_array(ghost_pa.get_carray('v'), shift, start=start) + if periodic_in_z and self.gamma_zx: + shift = xtranslate*self.gamma_zx + L = self.zmin + U = self.zmax + self._shift_periodic(ghost_pa.get_carray('z'), + shift*self.t, L, U, start=start) + self._add_to_array(ghost_pa.get_carray('w'), shift, start=start) # x_high start = ghost_pa.get_number_of_particles() @@ -870,29 +995,60 @@ cdef class CPUDomainManager(DomainManagerBase): ) self._add_to_array(ghost_pa.get_carray('x'), -xtranslate, start=start) + if periodic_in_y and self.gamma_yx: + shift = -xtranslate*self.gamma_yx + L = self.ymin + U = self.ymax + self._shift_periodic(ghost_pa.get_carray('y'), + shift*self.t, L, U, start=start) + self._add_to_array(ghost_pa.get_carray('v'), shift, start=start) + if periodic_in_z and self.gamma_zx: + shift = -xtranslate*self.gamma_zx + L = self.zmin + U = self.zmax + self._shift_periodic(ghost_pa.get_carray('z'), + shift*self.t, L, U, start=start) + self._add_to_array(ghost_pa.get_carray('w'), shift, start=start) if periodic_in_y: # Now do the corners from the previous. low.reset(); high.reset() - np = x.length + np = y.length for i in range(np): yi = y.data[i] if ( (yi - ymin) <= cell_size ): low.append(i) if ( (ymax - yi) <= cell_size ): high.append(i) + if np > 0: + start = ghost_pa.get_number_of_particles() + ghost_pa.extract_particles( + low, ghost_pa, align=False, props=copy_props[array_index] + ) + self._add_to_array(ghost_pa.get_carray('y'), ytranslate, + start=start) + if periodic_in_z and self.gamma_zy: + shift = ytranslate*self.gamma_zy + L = self.zmin + U = self.zmax + self._shift_periodic(ghost_pa.get_carray('z'), + shift*self.t, L, U, start=start) + self._add_to_array(ghost_pa.get_carray('w'), + shift, start=start) + + start = ghost_pa.get_number_of_particles() + ghost_pa.extract_particles( + high, ghost_pa, align=False, props=copy_props[array_index] + ) + self._add_to_array(ghost_pa.get_carray('y'), -ytranslate, + start=start) + if periodic_in_z and self.gamma_zy: + shift = ytranslate*self.gamma_zy + L = self.zmin + U = self.zmax + self._shift_periodic(ghost_pa.get_carray('z'), + -shift*self.t, L, U, start=start) + self._add_to_array(ghost_pa.get_carray('w'), + -shift, start=start) - start = ghost_pa.get_number_of_particles() - ghost_pa.extract_particles( - low, ghost_pa, align=False, props=copy_props[array_index] - ) - self._add_to_array(ghost_pa.get_carray('y'), ytranslate, - start=start) - - start = ghost_pa.get_number_of_particles() - ghost_pa.extract_particles( - high, ghost_pa, align=False, props=copy_props[array_index] - ) - self._add_to_array(ghost_pa.get_carray('y'), -ytranslate, - start=start) # Add the actual y_high and y_low now. # y_high @@ -902,6 +1058,15 @@ cdef class CPUDomainManager(DomainManagerBase): ) self._add_to_array(ghost_pa.get_carray('y'), -ytranslate, start=start) + if periodic_in_z and self.gamma_zy: + shift = ytranslate*self.gamma_zy + L = self.zmin + U = self.zmax + self._shift_periodic(ghost_pa.get_carray('z'), + -shift*self.t, L, U, start=start) + self._add_to_array(ghost_pa.get_carray('w'), + -shift, start=start) + # added.append_parray(copy) # y_low start = ghost_pa.get_number_of_particles() @@ -910,29 +1075,36 @@ cdef class CPUDomainManager(DomainManagerBase): ) self._add_to_array(ghost_pa.get_carray('y'), ytranslate, start=start) + if periodic_in_z and self.gamma_zy: + shift = ytranslate*self.gamma_zy + L = self.zmin + U = self.zmax + self._shift_periodic(ghost_pa.get_carray('z'), shift*self.t, L, U, start=start) + self._add_to_array(ghost_pa.get_carray('w'), shift, start=start) if periodic_in_z: # Now do the corners from the previous. low.reset(); high.reset() - np = x.length + np = z.length for i in range(np): zi = z.data[i] if ( (zi - zmin) <= cell_size ): low.append(i) if ( (zmax - zi) <= cell_size ): high.append(i) - start = ghost_pa.get_number_of_particles() - ghost_pa.extract_particles( - low, ghost_pa, align=False, props=copy_props[array_index] - ) - self._add_to_array(ghost_pa.get_carray('z'), ztranslate, - start=start) - - start = ghost_pa.get_number_of_particles() - ghost_pa.extract_particles( - high, ghost_pa, align=False, props=copy_props[array_index] - ) - self._add_to_array(ghost_pa.get_carray('z'), -ztranslate, - start=start) + if np > 0: + start = ghost_pa.get_number_of_particles() + ghost_pa.extract_particles( + low, ghost_pa, align=False, props=copy_props[array_index] + ) + self._add_to_array(ghost_pa.get_carray('z'), ztranslate, + start=start) + + start = ghost_pa.get_number_of_particles() + ghost_pa.extract_particles( + high, ghost_pa, align=False, props=copy_props[array_index] + ) + self._add_to_array(ghost_pa.get_carray('z'), -ztranslate, + start=start) # Add the actual z_high and z_low now. # z_high diff --git a/pysph/base/utils.py b/pysph/base/utils.py index 038f1b2be..393a6658d 100644 --- a/pysph/base/utils.py +++ b/pysph/base/utils.py @@ -526,3 +526,117 @@ def is_overloaded_method(method): break return count > 1 + + +def get_particle_array_beadchain_solid(constants=None, **props): + """Return a particle array for the BeadChainScheme for a fluid and solid. + + Parameters + ---------- + constants : dict + Dictionary of constants + + Other Parameters + ---------------- + props : dict + Additional keywords passed are set as the property arrays. + + See Also + -------- + get_particle_array + + """ + solid_props = ['V', 'vmag2', 'wij', + 'ug', 'vg', 'wg', + 'uf', 'vf', 'wf', + 'auhat', 'avhat', 'awhat', + 'uhat', 'vhat', 'what', + 'dudx', 'dudy', 'dudz', + 'dvdx', 'dvdy', 'dvdz', + 'dwdx', 'dwdy', 'dwdz'] + + pa = get_particle_array( + constants=constants, additional_props=solid_props, **props + ) + pa.set_output_arrays(['x', 'y', 'z', 'u', 'v', 'w', 'rho', 'm', 'h', 'p', + 'pid', 'gid', 'V', 'ug', 'vg', 'wg']) + + return pa + + +def get_particle_array_beadchain_fluid(constants=None, **props): + """Return a particle array for the BeadChainScheme for a fluid and solid. + + Parameters + ---------- + constants : dict + Dictionary of constants + + Other Parameters + ---------------- + props : dict + Additional keywords passed are set as the property arrays. + + See Also + -------- + get_particle_array + + """ + fluid_props = ['V', 'vmag2', + 'ug', 'vg', 'wg', + 'auhat', 'avhat', 'awhat', + 'uhat', 'vhat', 'what' + ] + + pa = get_particle_array( + constants=constants, additional_props=fluid_props, **props + ) + pa.set_output_arrays(['x', 'y', 'z', 'u', 'v', 'w', 'rho', 'm', 'h', 'p', + 'pid', 'gid', 'V', 'tag']) + + return pa + + +def get_particle_array_beadchain_fiber(constants=None, **props): + """Return a particle array for the BeadChainScheme for a fiber. + + Parameters + ---------- + constants : dict + Dictionary of constants + + Other Parameters + ---------------- + props : dict + Additional keywords passed are set as the property arrays. + + See Also + -------- + get_particle_array + + """ + fiber_props = ['V', 'wf', 'uf', 'vf', 'wg', 'wij', 'vg', 'ug', 'phifrac', + 'awhat', 'avhat', 'auhat', 'vhat', 'what', 'uhat', 'vmag2', + 'phi0', 'fractag', 'rho0', 'holdtag', 'eu', 'ev', + 'ew', 'dudx', 'dudy', 'dudz', 'dvdx', 'dvdy', 'dvdz', + 'dwdx', 'dwdy', 'dwdz', 'Fx', 'Fy', 'Fz', 'ex', + 'ey', 'ez', 'lprev', 'lnext', + 'rxnext', 'rynext', 'rznext', 'rnext', 'rxprev', 'ryprev', + 'rzprev', 'rprev', 'fidx'] + + pa = get_particle_array( + constants=constants, additional_props=fiber_props, **props + ) + pa.set_output_arrays(['x', 'y', 'z', + 'u', 'v', 'w', + 'fidx', 'tag', 'rho', 'm', 'h', 'p', + 'holdtag', 'fractag', 'gid', 'V', 'pid', + 'Fx', 'Fy', 'Fz', + 'rxnext', 'rynext', 'rznext', 'rnext', + 'rxprev', 'ryprev', 'rzprev', 'rprev', + 'lprev', 'lnext', + 'dudx', 'dudy', 'dudz', + 'dvdx', 'dvdy', 'dvdz', + 'dwdx', 'dwdy', 'dwdz']) + + return pa diff --git a/pysph/examples/couette_lees_edwards.py b/pysph/examples/couette_lees_edwards.py new file mode 100644 index 000000000..ad4a3850e --- /dev/null +++ b/pysph/examples/couette_lees_edwards.py @@ -0,0 +1,123 @@ +"""Couette flow with Lees-Edwards boundary conditions. + +This is using the transport velocity formulation and Lees-Edwards Boundary +Conditions (300 seconds) to form a shearflow without walls. The shearflow is +superposed with a constant lateral velocity to test particle transition across +the boundary. +""" + +import numpy as np + +# PySPH imports +from pysph.base.nnps import DomainManager +from pysph.base.utils import get_particle_array +from pysph.solver.application import Application +from pysph.sph.scheme import TVFScheme + +# domain and reference values +Re = 0.0125 +d = 0.5 +Lx = 2 * d +Ly = 0.4 * Lx +rho0 = 1.0 +nu = 0.01 + +# compute shear rate +Vmax = nu * Re / (2 * d) +gamma = 2 * Vmax / Lx + +# compute reference pressure +c0 = 10 * Vmax +p0 = c0 * c0 * rho0 + +# numerical setup +dx = 0.05 +hdx = 1.0 + +# adaptive time steps +h0 = hdx * dx +dt_cfl = 0.25 * h0 / (c0 + Vmax) +dt_viscous = 0.125 * h0 ** 2 / nu +dt_force = 1.0 + +# time integration +tf = 3200.0 +dt = min(dt_cfl, dt_viscous, dt_force) + + +class CouetteFlow(Application): + def create_domain(self): + """Create a DomainManager with Lees-Edwards BCs. + + Both directions are set as periodic and the value gamma_yx ist set to + the shear rate to shift particles crossing the x-boundary. As the BC + has to keep track of time, the time step is passed as well. + """ + return DomainManager( + xmin=0, + xmax=Lx, + periodic_in_x=True, + ymin=0, + ymax=Ly, + periodic_in_y=True, + gamma_yx=gamma, + n_layers=1, + dt=dt, + ) + + def create_particles(self): + _x = np.arange(dx / 2, Lx, dx) + + # create the fluid particles + _y = np.arange(dx / 2, Ly, dx) + + x, y = np.meshgrid(_x, _y) + fx = x.ravel() + fy = y.ravel() + + fluid = get_particle_array( + name="fluid", x=fx, y=fy, rho=rho0 * np.ones_like(fx) + ) + + print("Couette flow :: Re = %g, dt = %g" % (Re, dt)) + + self.scheme.setup_properties([fluid]) + + # setup the particle properties + volume = dx * dx + + # mass is set to get the reference density of rho0 + fluid.m[:] = volume * rho0 + + # volume is set as dx^2 + fluid.V[:] = 1.0 / volume + + # smoothing lengths + fluid.h[:] = hdx * dx + + # initial speed + fluid.v[:] = (fluid.x[:] - Lx / 2) * gamma + fluid.u[:] = Vmax + + # return the particle list + return [fluid] + + def create_scheme(self): + s = TVFScheme( + ["fluid"], + [], + dim=2, + rho0=rho0, + c0=c0, + nu=nu, + p0=p0, + pb=p0, + h0=dx * hdx + ) + s.configure_solver(tf=tf, dt=dt, output_only_real=False) + return s + + +if __name__ == "__main__": + app = CouetteFlow() + app.run() diff --git a/pysph/examples/fiber/beam_collision.py b/pysph/examples/fiber/beam_collision.py new file mode 100644 index 000000000..302d0a6c0 --- /dev/null +++ b/pysph/examples/fiber/beam_collision.py @@ -0,0 +1,243 @@ +"""Collision of a fiber in a damped field (10 minutes). + +################################################################################ +Beam Collision +################################################################################ + +Reference +--------- +N. Meyer et. al "Parameter Identification of Fiber Orientation Models Based on +Direct Fiber Simulation with Smoothed Particle Hydrodynamics", +Journal of Composites Science, 2020, 4, 77; doi:10.3390/jcs4020077 +""" +from math import sqrt + +import numpy as np + +from pysph.base.kernels import QuinticSpline +# PySPH imports +from pysph.base.utils import get_particle_array_beadchain_fiber +from pysph.solver.application import Application +from pysph.solver.solver import Solver +from pysph.sph.equation import Group +from pysph.sph.fiber.beadchain import Bending, Tension +from pysph.sph.fiber.utils import ComputeDistance, Contact, Damping, HoldPoints +from pysph.sph.integrator import EPECIntegrator +from pysph.sph.integrator_step import TransportVelocityStep +from pysph.sph.wc.transport_velocity import MomentumEquationPressureGradient + + +class Beam(Application): + def add_user_options(self, group): + group.add_argument( + "--D", action="store", type=float, dest="d", + default=1, help="Factor for damping. 1 is aperiodic limit." + ) + group.add_argument( + "--E", action="store", type=float, dest="E", + default=1E8, help="Young's modulus." + ) + group.add_argument( + "--N", action="store", type=int, dest="N", + default=10, help="Number of particles." + ) + group.add_argument( + "--gx", action="store", type=float, dest="gx", + default=0, help="Body force in x-direction." + ) + group.add_argument( + "--gy", action="store", type=float, dest="gy", + default=100, help="Body force in y-direction." + ) + group.add_argument( + "--gz", action="store", type=float, dest="gz", + default=0, help="Body force in z-direction." + ) + group.add_argument( + "--k", action="store", type=float, dest="k", + default=0.0, help="Friction coefficient." + ) + group.add_argument( + "--eta", action="store", type=float, dest="eta", + default=1.0, help="Absolute Viscosity." + ) + + def consume_user_options(self): + # fiber length + self.reducedL = 10.0 + + # numerical setup + self.N = self.options.N + self.reducedL = self.reducedL/(1-1/(2*self.N)) + self.dx = self.reducedL/self.N # particle spacing + self.h = self.dx + + # fluid properties + self.rho0 = 1.0 + self.p0 = 1.0 + + # fiber properties + self.A = 1.0 + self.Ip = self.A/12.0 + self.E = self.options.E + # Analytical solution for angular eigenfrequencies: + # Pi/L np.sqrt(E/rho) (2n-1)/2 + # --> first analytical eigenfrequency: + self.omega0_tension = np.pi/(2*self.reducedL)*np.sqrt(self.E/self.rho0) + self.omega0_bending = 3.5156*np.sqrt( + self.E*self.Ip/(self.rho0*self.A*self.reducedL**4)) + if self.options.gx > self.options.gy: + self.omega0 = self.omega0_tension + else: + self.omega0 = self.omega0_bending + m = self.rho0*self.A*self.dx + self.D = self.options.d*m*self.omega0 + self.gx = self.options.gx + self.gy = self.options.gy + self.gz = self.options.gz + print('Damping: %g, Omega0: %g' % (self.D, self.omega0)) + + # setup time step + if abs(self.gx) > 0.0 or abs(self.gy) > 0.0 or abs(self.gz) > 0.0: + dt_force = 0.25 * np.sqrt( + self.h/(sqrt(self.gx**2+self.gy**2+self.gz**2))) + else: + dt_force = 10000 + dt_tension = 0.5*self.h*np.sqrt(self.rho0/self.E) + dt_bending = 0.5*self.h**2*np.sqrt(self.rho0*self.A/(self.E*2*self.Ip)) + + self.tf = 20 + + self.dt = min(dt_force, dt_tension, dt_bending) + + def create_scheme(self): + return None + + def create_particles(self): + _x = np.linspace(-self.dx, self.reducedL-self.dx, self.N+1) + _y = np.array([0.0]) + _z = np.array([0.0]) + x, y, z = np.meshgrid(_x, _y, _z) + fiber1_x = x.ravel() + fiber1_y = y.ravel() + fiber1_z = z.ravel() + + _x = np.array([0.75*self.reducedL]) + _y = np.array([-2*self.dx]) + _z = -np.linspace(-0.25*self.reducedL, 0.75*self.reducedL, self.N+1) + # _x = np.linspace(0.0, self.reducedL, self.N+1) + # _y = np.array([-1.5*self.dx]) + # _z = np.array([0.0]) + x, y, z = np.meshgrid(_x, _y, _z) + fiber2_x = x.ravel() + fiber2_y = y.ravel() + fiber2_z = z.ravel() + + # volume is set as dx * A + volume = self.A * self.dx + + # create arrays + fiber1 = get_particle_array_beadchain_fiber( + name='fiber1', x=fiber1_x, y=fiber1_y, z=fiber1_z, + m=volume*self.rho0, rho=self.rho0, h=self.h, lprev=self.dx, + lnext=self.dx, phi0=np.pi, phifrac=2.0, fidx=range(self.N+1), + V=1./volume) + fiber2 = get_particle_array_beadchain_fiber( + name='fiber2', x=fiber2_x, y=fiber2_y, z=fiber2_z, + m=volume*self.rho0, rho=self.rho0, h=self.h, lprev=self.dx, + lnext=self.dx, phi0=np.pi, phifrac=2.0, fidx=range(self.N+1), + V=1./volume) + + # tag particles to be hold + fiber1.holdtag[:] = 0 + fiber1.holdtag[0] = 2 + fiber1.holdtag[1] = 1 + fiber1.holdtag[2] = 2 + + # return the particle list + return [fiber1, fiber2] + + def create_equations(self): + equations = [ + Group( + equations=[ + ComputeDistance(dest='fiber1', sources=['fiber1']), + ComputeDistance(dest='fiber2', sources=['fiber2']), + ], + ), + Group( + equations=[ + MomentumEquationPressureGradient( + dest='fiber1', + sources=['fiber1', 'fiber2'], pb=0.0, + gx=self.gx, gy=self.gy, gz=self.gz), + MomentumEquationPressureGradient( + dest='fiber2', + sources=['fiber1', 'fiber2'], pb=0.0, + gx=self.gx, gy=self.gy, gz=self.gz), + Tension( + dest='fiber1', + sources=None, + ea=self.E*self.A), + Tension( + dest='fiber2', + sources=None, + ea=self.E*self.A), + Bending( + dest='fiber1', + sources=None, + ei=self.E*self.Ip), + Bending( + dest='fiber2', + sources=None, + ei=self.E*self.Ip), + Contact( + dest='fiber1', + sources=['fiber1', 'fiber2'], + E=self.E, d=self.dx, dim=3, + k=self.options.k, + eta0=self.options.eta, + dt=self.dt), + Contact( + dest='fiber2', + sources=['fiber1', 'fiber2'], + E=self.E, d=self.dx, dim=3, + k=self.options.k, + eta0=self.options.eta, + dt=self.dt), + Damping( + dest='fiber1', + sources=None, + d=self.D), + Damping( + dest='fiber2', + sources=None, + d=self.D) + ], + ), + Group( + equations=[ + HoldPoints(dest='fiber1', sources=None, tag=2, x=False), + HoldPoints(dest='fiber1', sources=None, tag=1, y=False), + HoldPoints(dest='fiber2', sources=None, tag=2, x=False), + HoldPoints(dest='fiber2', sources=None, tag=1, y=False), + ], + ), + ] + return equations + + def create_solver(self): + """Set up the default integrator for fiber particles.""" + kernel = QuinticSpline(dim=3) + integrator = EPECIntegrator( + fiber1=TransportVelocityStep(), + fiber2=TransportVelocityStep()) + solver = Solver( + kernel=kernel, dim=3, integrator=integrator, dt=self.dt, + tf=self.tf, N=200) + return solver + + +if __name__ == '__main__': + app = Beam() + app.run() diff --git a/pysph/examples/fiber/rve.py b/pysph/examples/fiber/rve.py new file mode 100644 index 000000000..7cddcfb03 --- /dev/null +++ b/pysph/examples/fiber/rve.py @@ -0,0 +1,483 @@ +"""Example for Mini RVE. + +################################################################################ +Mini RVE +################################################################################ + +Reference +--------- +N. Meyer et. al "Parameter Identification of Fiber Orientation Models Based on +Direct Fiber Simulation with Smoothed Particle Hydrodynamics", +Journal of Composites Science, 2020, 4, 77; doi:10.3390/jcs4020077 +""" +import itertools +# general imports +import os +import random + +import numpy as np + +from pysph.base.nnps import DomainManager +from pysph.base.utils import (get_particle_array_beadchain_fiber, + get_particle_array_beadchain_fluid) +from pysph.solver.application import Application +from pysph.solver.utils import load, remove_irrelevant_files +from pysph.sph.scheme import BeadChainScheme + +# from pysph.solver.tools import FiberIntegrator + + +class RVE(Application): + """Generate a mini RVE and evaluate its fiber orientation tensor.""" + + def create_scheme(self): + """Use the BeadChainScheme for this model.""" + return BeadChainScheme(["fluid"], [], ["fibers"], dim=3) + + def add_user_options(self, group): + """Set these parameters in the command line.""" + group.add_argument( + "--dx", + action="store", + type=float, + dest="dx", + default=0.0001, + help="Particle spacing", + ) + group.add_argument( + "--lf", + action="store", + type=int, + dest="lf", + default=5, + help="Fiber length in multiples of dx", + ) + group.add_argument( + "--mu", + action="store", + type=float, + dest="mu", + default=1.0, + help="Absolute viscosity", + ) + group.add_argument( + "--S", + action="store", + type=float, + dest="S", + default=100, + help="Dimensionless fiber stiffness", + ) + group.add_argument( + "--G", + action="store", + type=float, + dest="G", + default=1.0, + help="Shear rate applied to the cube", + ) + group.add_argument( + "--Re", + action="store", + type=float, + dest="Re", + default=0.5, + help="Desired Reynolds number.", + ) + group.add_argument( + "--volfrac", + action="store", + type=float, + dest="vol_frac", + default=0.0014, + help="Volume fraction of fibers in suspension.", + ) + group.add_argument( + "--rot", + action="store", + type=float, + dest="rot", + default=2.0, + help="Number of half rotations.", + ) + group.add_argument( + "--C", + action="store", + type=float, + dest="C", + default=15.0, + help="Cube size as multiples of particle spacing.", + ) + group.add_argument( + "--continue", + action="store", + type=str, + dest="continuation", + default=None, + help="Set a file for continuation of run.", + ) + + def consume_user_options(self): + """Initialize geometry, properties and time stepping.""" + # Initial spacing of particles + self.dx = self.options.dx + self.h0 = self.dx + + # The fiber length is the aspect ratio times fiber diameter + self.L = self.options.lf * self.dx + + # Cube size + self.C = self.options.C * self.dx + + # Density from Reynolds number + self.Vmax = self.options.G * self.C / 2.0 + self.rho0 = (self.options.mu * self.options.Re) / (self.Vmax * self.C) + + # The kinematic viscosity + self.nu = self.options.mu / self.rho0 + + # Fiber aspect ratio (assuming cylindrical shape) + R = self.dx / (np.sqrt(np.pi)) + self.d = 2.0 * R + self.ar = self.L / self.d + print("Aspect ratio is %f" % self.ar) + + # cross section properties + self.A = np.pi * R ** 2.0 + self.Ip = np.pi * R ** 4.0 / 4.0 + + # inertia properties + mass = 3.0 * self.rho0 * self.dx * self.A + self.J = ( + 1.0 / 4.0 * mass * R ** 2.0 + + 1.0 / 12.0 * mass * (3.0 * self.dx) ** 2.0 + ) + + # stiffness from dimensionless stiffness + self.E = ( + 4.0 / np.pi * ( + self.options.S * self.options.mu * self.options.G * self.ar) + ) + + # The speed of sound c0 is computed as 10 times the maximum velocity. + # This should keep the density change within 1% + self.c0 = 10.0 * self.Vmax + self.p0 = self.c0 ** 2 * self.rho0 + + # Background pressure in Adami's transport velocity formulation + self.pb = self.p0 + + # Simulation time + self.t = ( + self.options.rot * np.pi * (self.ar + 1.0 / self.ar) + / self.options.G) + print("Simulated time is %g s" % self.t) + + # Determine number of fibers to be generated + vol_fiber = self.L * self.dx * self.dx + vol = self.C ** 3 + self.n = int(round(self.options.vol_frac * vol / vol_fiber)) + + def configure_scheme(self): + """Set up scheme and solver. + + The flag 'direct' means that elastic equations and contact equations + are solved together with all other equations. If the fiber is very + stiff, one may use a subcycle to integrate fiber positions. Therefore, + set 'direct=False' and uncomment the FiberIntegrator tool. + """ + self.scheme.configure( + rho0=self.rho0, + c0=self.c0, + nu=self.nu, + p0=self.p0, + pb=self.pb, + h0=self.h0, + dx=self.dx, + A=self.A, + Ip=self.Ip, + J=self.J, + E=self.E, + d=self.d, + direct=True, + ) + if self.n < 1: + self.scheme.configure(fibers=[]) + self.scheme.configure_solver( + tf=self.t, + # pfreq=1, + N=self.options.rot * 100, + ) + + def create_particles(self): + """Create or load particle arrays.""" + if self.options.continuation: + data = load(self.options.continuation) + fluid = data["arrays"]["fluid"] + fibers = data["arrays"]["fibers"] + fibers.phifrac[:] = 2.0 + fibers.phi0[:] = np.pi + self.solver.t = data["solver_data"]["t"] + self.solver.count = data["solver_data"]["count"] + return [fluid, fibers] + else: + return self.create_suspension_particles() + + def create_suspension_particles(self): + """Create particle arrays.""" + fdx = self.dx + dx2 = fdx / 2 + + # Computation of each particles initial volume. + volume = fdx ** 3 + + # Mass is set to get the reference density of rho0. + mass = volume * self.rho0 + + # Initial inverse volume (necessary for transport velocity equations) + V = 1.0 / volume + + # Create grid points for particles + _x = np.arange(dx2, self.C, fdx) + _y = np.arange(dx2, self.C, fdx) + _z = np.arange(dx2, self.C, fdx) + fx, fy, fz = self.get_meshgrid(_x, _y, _z) + + # Remove particles at fiber position and create fiber particles. + indices = [] + fibers = [] + fibx = tuple() + fiby = tuple() + fibz = tuple() + + positions = list(itertools.product(_x, _y, _z)) + random.shuffle(positions) + N = 0 + while N < self.n: + xx, yy, zz = positions.pop() + idx_list = [] + for i in range(len(fx)): + # periodic extending above + if xx + self.L / 2 > self.C: + if ( + (fx[i] < (xx + self.L / 2 - self.C) + or fx[i] > xx - self.L / 2) + and fy[i] < yy + self.dx / 2 + and fy[i] > yy - self.dx / 2 + and fz[i] < zz + self.dx / 2 + and fz[i] > zz - self.dx / 2 + ): + idx_list.append(i) + # periodic extending below + elif xx - self.L / 2 < 0: + if ( + (fx[i] < xx + self.L / 2 + or fx[i] > (xx - self.L / 2 + self.C)) + and fy[i] < yy + self.dx / 2 + and fy[i] > yy - self.dx / 2 + and fz[i] < zz + self.dx / 2 + and fz[i] > zz - self.dx / 2 + ): + idx_list.append(i) + # standard case + else: + if ( + fx[i] < xx + self.L / 2 + and fx[i] > xx - self.L / 2 + and fy[i] < yy + self.dx / 2 + and fy[i] > yy - self.dx / 2 + and fz[i] < zz + self.dx / 2 + and fz[i] > zz - self.dx / 2 + ): + idx_list.append(i) + + idx_set = set(idx_list) + if len(idx_set.intersection(set(indices))) == 0: + N = N + 1 + indices = indices + idx_list + + # Generate fiber particles + if self.options.lf % 2 == 1: + _fibx = np.linspace( + xx - self.options.lf // 2 * self.dx, + xx + self.options.lf // 2 * self.dx, + self.options.lf, + ) + else: + _fibx = np.arange( + xx - self.L / 2, xx + self.L / 2 - self.dx / 4, self.dx + ) + _fiby = np.array([yy]) + _fibz = np.array([zz]) + _fibx, _fiby, _fibz = self.get_meshgrid(_fibx, _fiby, _fibz) + fibx = fibx + (_fibx,) + fiby = fiby + (_fiby,) + fibz = fibz + (_fibz,) + + print("Created %d fibers." % N) + + # Finally create all particle arrays. Note that fluid particles are + # removed in the area, where the fiber is placed. + fluid = get_particle_array_beadchain_fluid( + name="fluid", + x=fx, + y=fy, + z=fz, + m=mass, + rho=self.rho0, + h=self.h0, + V=V + ) + fluid.remove_particles(indices) + + if self.n > 0: + fibers = get_particle_array_beadchain_fiber( + name="fibers", + x=np.concatenate(fibx), + y=np.concatenate(fiby), + z=np.concatenate(fibz), + m=mass, + rho=self.rho0, + h=self.h0, + lprev=self.dx, + lnext=self.dx, + phi0=np.pi, + phifrac=2.0, + fidx=range(self.options.lf * self.n), + V=V, + ) + # 'Break' fibers in segments + endpoints = [i * self.options.lf - 1 for i in range(1, self.n)] + fibers.fractag[endpoints] = 1 + + # Setting the initial velocities for a shear flow. + fluid.v[:] = self.options.G * (fluid.x[:] - self.C / 2) + + fibers.v[:] = self.options.G * (fibers.x[:] - self.C / 2) + return [fluid, fibers] + + def create_domain(self): + """Create periodic boundary conditions in all directions. + + Additionally, gamma values are set to enforce Lee-Edwards BCs. + """ + return DomainManager( + xmin=0, + xmax=self.C, + periodic_in_x=True, + ymin=0, + ymax=self.C, + periodic_in_y=True, + zmin=0, + zmax=self.C, + periodic_in_z=True, + gamma_yx=self.options.G, + n_layers=1, + dt=self.solver.dt, + calls_per_step=2, + ) + + # def create_tools(self): + # """Set up fiber integrator.""" + # return [FiberIntegrator(self.particles, self.scheme, self.domain, + # D=0.002*self.options.lf)] + + def get_meshgrid(self, xx, yy, zz): + """Create meshgrids.""" + x, y, z = np.meshgrid(xx, yy, zz) + x = x.ravel() + y = y.ravel() + z = z.ravel() + return [x, y, z] + + def post_process(self, info_fname): + """Save fiber orientation tensor to csv file.""" + if len(self.output_files) == 0: + return + + from pysph.tools.pprocess import get_ke_history + from matplotlib import pyplot as plt + + t, ke = get_ke_history(self.output_files, "fluid") + plt.clf() + plt.plot(t, ke) + plt.xlabel("t") + plt.ylabel("Kinetic energy") + fig = os.path.join(self.output_dir, "ke_history.png") + plt.savefig(fig, dpi=300) + + # empty list for time + t = [] + + # empty lists for fiber orientation tensors + A = [] + + # iteration over all output files + output_files = remove_irrelevant_files(self.output_files) + for fname in output_files: + data = load(fname) + + # extracting time + t.append(data["solver_data"]["t"]) + + if self.n > 0: + # extrating all arrays. + directions = [] + fiber = data["arrays"]["fibers"] + startpoints = [ + i * (self.options.lf - 1) for i in range(0, self.n) + ] + endpoints = [ + i * (self.options.lf - 1) - 1 for i in range(1, self.n + 1) + ] + for start, end in zip(startpoints, endpoints): + px = np.mean(fiber.rxnext[start:end]) + py = np.mean(fiber.rynext[start:end]) + pz = np.mean(fiber.rznext[start:end]) + + n = np.array([px, py, pz]) + norm = np.linalg.norm(n) + if norm == 0: + p = np.array([1, 0, 0]) + else: + p = n / norm + directions.append(p) + + N = len(directions) + a = np.zeros([3, 3]) + for p in directions: + for i in range(3): + for j in range(3): + a[i, j] += 1.0 / N * (p[i] * p[j]) + A.append(a.ravel()) + + csv_file = os.path.join(self.output_dir, "A.csv") + data = np.hstack((np.matrix(t).T, np.vstack(A))) + np.savetxt(csv_file, data, delimiter=",") + + # plot results + data = np.loadtxt(csv_file, delimiter=',') + t = data[:, 0] + plt.figure(figsize=(12, 3)) + for j, i in enumerate([0, 4, 8, 1]): + plt.subplot("14"+str(j+1)) + p = plt.plot(self.options.G*t, data[:, i+1]) + plt.xlabel("Strains") + plt.ylabel("Component %d" % j) + if i % 2 == 0: + plt.ylim([0, 1]) + else: + plt.ylim([-1, 1]) + plt.tight_layout() + plt.savefig(csv_file.replace('.csv', '.pdf')) + try: + from matplotlib2tikz import save as tikz_save + tikz_save(csv_file.replace('.csv', '.tex')) + except ImportError: + print("Did not write tikz figure.") + + +if __name__ == "__main__": + app = RVE() + app.run() + app.post_process(app.info_filename) diff --git a/pysph/examples/fiber/shearflow.py b/pysph/examples/fiber/shearflow.py new file mode 100644 index 000000000..38dc45754 --- /dev/null +++ b/pysph/examples/fiber/shearflow.py @@ -0,0 +1,361 @@ +"""Example for fiber shearflow. + +################################################################################ +3D shearflow with a single fiber +################################################################################ + +Reference +--------- +N. Meyer et. al "Parameter Identification of Fiber Orientation Models Based on +Direct fiber Simulation with Smoothed Particle Hydrodynamics", +Journal of Composites Science, 2020, 4, 77; doi:10.3390/jcs4020077 +""" + +import os + +import numpy as np + +from pysph.base.nnps import DomainManager +from pysph.base.utils import (get_particle_array_beadchain_fiber, + get_particle_array_beadchain_fluid, + get_particle_array_beadchain_solid) +from pysph.solver.application import Application +from pysph.solver.utils import load, remove_irrelevant_files +from pysph.sph.scheme import BeadChainScheme + +# from pysph.solver.tools import FiberIntegrator + + +def get_zhang_aspect_ratio(aspect_ratio): + """Jeffery's equivalent aspect ratio. + + Approximation for small aspect ratios from + Zhang et al. 2011 + """ + return (0.000035*aspect_ratio**3 - 0.00467*aspect_ratio**2 + + 0.764*aspect_ratio + 0.404) + + +def jeffery_ode(phi, t, ar, G): + """Jeffery's Equation for planar rotation of a rigid.""" + lbd = (ar**2. - 1.)/(ar**2. + 1.) + return 0.5*G*(1. + lbd*np.cos(2.*phi)) + + +class Channel(Application): + """Application for the channel flow driven by top an bottom walls.""" + + def create_scheme(self): + """Use BeadChainScheme for this application.""" + return BeadChainScheme(['fluid'], ['channel'], ['fiber'], dim=3) + + def add_user_options(self, group): + """Add options to aplication.""" + group.add_argument( + "--dx", action="store", type=float, dest="dx", + default=0.0001, help="Particle Spacing" + ) + group.add_argument( + "--lf", action="store", type=int, dest="lf", + default=5, help="Fiber length in multiples of dx" + ) + group.add_argument( + "--mu", action="store", type=float, dest="mu", + default=1.0, help="Absolute viscosity" + ) + group.add_argument( + "--S", action="store", type=float, dest="S", + default=100, help="Dimensionless fiber stiffness" + ) + group.add_argument( + "--G", action="store", type=float, dest="G", + default=1.0, help="Shear rate" + ) + group.add_argument( + "--Re", action="store", type=float, dest="Re", + default=0.5, help="Desired Reynolds number." + ) + group.add_argument( + "--frac", action="store", type=float, dest="phifrac", + default=5.0, help="Critical bending angle for fracture." + ) + group.add_argument( + "--rot", action="store", type=float, dest="rot", + default=1.0, help="Number of half rotations." + ) + + def consume_user_options(self): + """Initialize geometry, properties and time stepping.""" + # Initial spacing of particles + self.dx = self.options.dx + self.h0 = self.dx + + # The fiber length is the aspect ratio times fiber diameter + self.Lf = self.options.lf*self.dx + + # Use fiber aspect ratio to determine the channel width. + self.Ly = self.Lf + 2.*self.dx + + # Density from Reynolds number + self.Vmax = self.options.G*self.Ly/2. + self.rho0 = (self.options.mu*self.options.Re)/(self.Vmax*self.Lf) + + # The channel length is twice the width + dx to make it symmetric. + self.Lx = 2.*self.Ly + self.dx + + # The position of the fiber's center is set to the center of the + # channel. + self.x_fiber = 0.5*self.Lx + self.y_fiber = 0.5*self.Ly + self.z_fiber = 0.5*self.Ly + + # The kinematic viscosity is computed from absolute viscosity and + # scaled density. + self.nu = self.options.mu/self.rho0 + + # mass properties + R = self.dx/np.sqrt(np.pi) # Assuming cylindrical shape + self.d = 2.*R + self.ar = self.Lf/self.d # Actual fiber aspect ratio + print('Aspect ratio is %f' % self.ar) + + self.A = np.pi*R**2. + self.Ip = np.pi*R**4./4. + mass = 3.*self.rho0*self.dx*self.A + self.J = 1./4.*mass*R**2. + 1./12.*mass*(3.*self.dx)**2. + + # stiffness from dimensionless stiffness + self.E = 4.0/np.pi*( + self.options.S*self.options.mu*self.options.G*self.ar) + + # SPH uses weakly compressible fluids. Therefore, the speed of sound c0 + # is computed as 10 times the maximum velocity. This should keep the + # density change within 1% + self.c0 = 10.*self.Vmax + self.p0 = self.c0**2*self.rho0 + + # Background pressure in Adami's transport velocity formulation + self.pb = self.p0 + + # Time + ar = get_zhang_aspect_ratio(self.ar) + self.t = self.options.rot*np.pi*(ar + 1./ar)/self.options.G + print("Simulated time is %g s" % self.t) + + def configure_scheme(self): + """Set up solver and scheme.""" + self.scheme.configure( + rho0=self.rho0, + c0=self.c0, + nu=self.nu, + p0=self.p0, + pb=self.pb, + h0=self.h0, + dx=self.dx, + A=self.A, + Ip=self.Ip, + J=self.J, + E=self.E, + d=self.d, + direct=True) + self.scheme.configure_solver( + tf=self.t, + # pfreq=1, + N=self.options.rot*200 + ) + + def create_particles(self): + """Three particle arrays are created. + + A fluid, representing the polymer matrix, a fiber with additional + properties and a channel of dummyparticles. + """ + # The fluid might be scaled compared to the fiber. fdx is a shorthand + # for the fluid spacing and dx2 is a shorthand for the half of it. + fdx = self.dx + dx2 = fdx/2. + + # Creating grid points for particles + _x = np.arange(dx2, self.Lx, fdx) + _y = np.arange(dx2, self.Ly, fdx) + _z = np.arange(dx2, self.Ly, fdx) + fx, fy, fz = self.get_meshgrid(_x, _y, _z) + + # Remove particles at fiber position + indices = [] + for i in range(len(fx)): + xx = self.x_fiber + yy = self.y_fiber + zz = self.z_fiber + + # vertical + if (fx[i] < xx + self.dx/2. and fx[i] > xx - self.dx/2. and + fy[i] < yy + self.Lf/2. and fy[i] > yy - self.Lf/2. and + fz[i] < zz + self.dx/2. and fz[i] > zz - self.dx/2.): + indices.append(i) + + # create vertical fiber + _fibx = np.array([xx]) + _fiby = np.arange(yy - self.Lf/2. + self.dx/2., + yy + self.Lf/2. + self.dx/4., + self.dx) + + _fibz = np.array([zz]) + fibx, fiby, fibz = self.get_meshgrid(_fibx, _fiby, _fibz) + + # Determine the size of dummy region + ghost_extent = 3.*fdx + + # Create the channel particles at the top + _y = np.arange(self.Ly + dx2, self.Ly + dx2 + ghost_extent, fdx) + tx, ty, tz = self.get_meshgrid(_x, _y, _z) + + # Create the channel particles at the bottom + _y = np.arange(-dx2, -dx2 - ghost_extent, -fdx) + bx, by, bz = self.get_meshgrid(_x, _y, _z) + + # Concatenate the top and bottom arrays + cx = np.concatenate((tx, bx)) + cy = np.concatenate((ty, by)) + cz = np.concatenate((tz, bz)) + + # Computation of each particles initial volume. + volume = fdx**3. + + # Mass is set to get the reference density of rho0. + mass = volume*self.rho0 + + # assign unique ID (within fiber) to each fiber particle. + fidx = range(0, self.options.lf) + + # Initial inverse volume (necessary for transport velocity equations) + V = 1./volume + + # Finally create all particle arrays. Note that fluid particles are + # removed in the area, where the fiber is placed. + channel = get_particle_array_beadchain_solid( + name='channel', x=cx, y=cy, z=cz, m=mass, rho=self.rho0, + h=self.h0, V=V) + fluid = get_particle_array_beadchain_fluid( + name='fluid', x=fx, y=fy, z=fz, m=mass, rho=self.rho0, + h=self.h0, V=V) + fluid.remove_particles(indices) + fiber = get_particle_array_beadchain_fiber( + name='fiber', x=fibx, y=fiby, z=fibz, m=mass, + rho=self.rho0, h=self.h0, lprev=self.dx, lnext=self.dx, + phi0=np.pi, phifrac=self.options.phifrac, fidx=fidx, V=V) + + # The number of fiber particles should match the aspect ratio. This + # assertation fails, if something was wrong in the fiber generation. + assert(fiber.get_number_of_particles() == self.options.lf) + + # Setting the initial velocities for a shear flow. + fluid.u[:] = self.options.G*(fluid.y[:] - self.Ly/2.) + fiber.u[:] = self.options.G*(fiber.y[:] - self.Ly/2.) + channel.u[:] = self.options.G*np.sign(channel.y[:])*self.Ly/2. + + # Return the particle list. + return [fluid, channel, fiber] + + def create_domain(self): + """Create periodic boundary conditions in x-direction.""" + return DomainManager(xmin=0, xmax=self.Lx, zmin=0, zmax=self.Ly, + periodic_in_x=True, periodic_in_z=True) + + # def create_tools(self): + # """Add an integrator for the fiber.""" + # return [FiberIntegrator(self.particles, self.scheme, self.domain)] + + def get_meshgrid(self, xx, yy, zz): + """Generate meshgrids quickly.""" + x, y, z = np.meshgrid(xx, yy, zz) + x = x.ravel() + y = y.ravel() + z = z.ravel() + return [x, y, z] + + def _plots(self): + """Create plots. + + It is employing a iteration over all time steps. + """ + try: + from matplotlib import pyplot as plt + from scipy.integrate import odeint + except ImportError: + print("Postprocessing requires scipy and matplotlib.") + return + # empty list for time and orientation angle + t = [] + angle = [] + N = 0 + + # iteration over all output files + output_files = remove_irrelevant_files(self.output_files) + for i, fname in enumerate(output_files): + data = load(fname) + # extracting time and fiber data + t.append(data['solver_data']['t']) + fiber = data['arrays']['fiber'] + + # computation of orientation angle + dxx = fiber.x[0] - fiber.x[-1] + dyy = fiber.y[0] - fiber.y[-1] + a = np.arctan(dxx / (dyy + 0.01 * self.h0)) + N * np.pi + if len(angle) > 0 and a - angle[-1] > 3: + N -= 1 + a -= np.pi + elif len(angle) > 0 and a - angle[-1] < -3: + N += 1 + a += np.pi + angle.append(a) + + # Integrate Jeffery's solution + print("Solving Jeffery's ODE") + t = np.array(t) + phi0 = angle[0] + ar_zhang = get_zhang_aspect_ratio(self.ar) + angle_jeffery_zhang = odeint(jeffery_ode, phi0, t, atol=1E-15, + args=(ar_zhang, self.options.G)) + + # open new plot + plt.figure() + + # plot computed angle and Jeffery's solution + plt.plot(t*self.options.G, angle, '-k') + plt.plot(t*self.options.G, angle_jeffery_zhang, '--k', color='grey') + + # labels + plt.xlabel('Strains $tG$') + plt.ylabel(r'Rotation angle $\phi$') + plt.legend(['SPH Simulation', 'Jeffery (Zhang)']) + plt.grid() + x1, x2, y1, y2 = plt.axis() + plt.axis((0, x2, 0, y2)) + ax = plt.gca() + ax.set_yticks([0, 0.5*np.pi, np.pi, 1.5*np.pi]) + ax.set_yticklabels(['0', r'$\pi/2$', r'$\pi$', r'$3/2\pi$']) + plt.tight_layout() + + # save figure + angfig = os.path.join(self.output_dir, 'angleplot.pdf') + plt.savefig(angfig, dpi=300, bbox_inches='tight') + try: + tex_fig = os.path.join(self.output_dir, "angleplot.tex") + from tikzplotlib import save as tikz_save + tikz_save(tex_fig) + except ImportError: + print("Did not write tikz figure.") + print("Angleplot written to %s." % angfig) + + def post_process(self, info_fname): + """Build plots and files as results.""" + if len(self.output_files) == 0: + return + self._plots() + + +if __name__ == '__main__': + app = Channel() + app.run() + app.post_process(app.info_filename) diff --git a/pysph/solver/solver.py b/pysph/solver/solver.py index 3a32325f9..5798761dd 100644 --- a/pysph/solver/solver.py +++ b/pysph/solver/solver.py @@ -24,7 +24,8 @@ def __init__(self, dim=2, integrator=None, kernel=None, n_damp=0, tf=1.0, dt=1e-3, adaptive_timestep=False, cfl=0.3, output_at_times=(), - fixed_h=False, **kwargs): + fixed_h=False, + **kwargs): """**Constructor** Any additional keyword args are used to set the values of any @@ -113,6 +114,9 @@ def __init__(self, dim=2, integrator=None, kernel=None, # default output printing frequency self.pfreq = 100 + # default value for total number of dumps --> overrides pfreq + self.N = 0 + # Compress generated files. self.compress_output = False self.disable_output = False @@ -696,6 +700,8 @@ def _dump_output_if_needed(self): # dump output if the iteration number is a multiple of the printing # frequency. + if self.N > 0: + self.pfreq = int(self.tf/(self.N*self.dt)) dump = self.count % self.pfreq == 0 # Consider the other cases if user has requested output at a specified diff --git a/pysph/solver/tools.py b/pysph/solver/tools.py index 58f0af63b..174b55598 100644 --- a/pysph/solver/tools.py +++ b/pysph/solver/tools.py @@ -94,6 +94,174 @@ def post_step(self, solver): self.interp.nnps.update_domain() +class FiberIntegrator(Tool): + def __init__(self, all_particles, scheme, domain=None, innerloop=True, + updates=True, parallel=False, steps=None, D=0): + """The second integrator is a simple Euler-Integrator (accurate + enough due to very small time steps; very fast) using EBGSteps. + EBGSteps are basically the same as EulerSteps, exept for the fact + that they work with an intermediate ebg velocity [eu, ev, ew]. + This velocity does not interfere with the actual velocity, which + is neseccery to not disturb the real velocity through artificial + damping in this step. The ebg velocity is initialized for each + inner loop again and reset in the outer loop.""" + from math import ceil + from pysph.base.kernels import CubicSpline + from pysph.sph.integrator_step import EBGStep + from compyle.config import get_config + from pysph.sph.integrator import EulerIntegrator + from pysph.sph.scheme import BeadChainScheme + from pysph.sph.equation import Group + from pysph.sph.fiber.utils import (HoldPoints, Contact, + ComputeDistance) + from pysph.sph.fiber.beadchain import (Tension, Bending, + ArtificialDamping) + from pysph.base.nnps import DomainManager, LinkedListNNPS + from pysph.sph.acceleration_eval import AccelerationEval + from pysph.sph.sph_compiler import SPHCompiler + + if not isinstance(scheme, BeadChainScheme): + raise TypeError("Scheme must be BeadChainScheme") + + self.innerloop = innerloop + self.dt = scheme.dt + self.fiber_dt = scheme.fiber_dt + self.domain_updates = updates + self.steps = steps + self.D = D + self.eta0 = scheme.rho0 * scheme.nu + + # if there are more than 1 particles involved, elastic equations are + # iterated in an inner loop. + if self.innerloop: + # second integrator + # self.fiber_integrator = EulerIntegrator(fiber=EBGStep()) + steppers = {} + for f in scheme.fibers: + steppers[f] = EBGStep() + self.fiber_integrator = EulerIntegrator(**steppers) + # The type of spline has no influence here. It must be large enough + # to contain the next particle though. + kernel = CubicSpline(dim=scheme.dim) + equations = [] + g1 = [] + for fiber in scheme.fibers: + g1.append(ComputeDistance(dest=fiber, sources=[fiber])) + equations.append(Group(equations=g1)) + + g2 = [] + for fiber in scheme.fibers: + g2.append(Tension(dest=fiber, + sources=None, + ea=scheme.E*scheme.A)) + g2.append(Bending(dest=fiber, + sources=None, + ei=scheme.E*scheme.Ip)) + g2.append(Contact(dest=fiber, + sources=scheme.fibers, + E=scheme.E, + d=scheme.dx, + dim=scheme.dim, + k=scheme.k, + lim=scheme.lim, + eta0=self.eta0)) + g2.append(ArtificialDamping(dest=fiber, + sources=None, + d=self.D)) + equations.append(Group(equations=g2)) + + g3 = [] + for fiber in scheme.fibers: + g3.append(HoldPoints(dest=fiber, sources=None, tag=100)) + equations.append(Group(equations=g3)) + + # These equations are applied to fiber particles only - that's the + # reason for computational speed up. + particles = [p for p in all_particles if p.name in scheme.fibers] + # A seperate DomainManager is needed to ensure that particles don't + # leave the domain. + if domain: + xmin = domain.manager.xmin + ymin = domain.manager.ymin + zmin = domain.manager.zmin + xmax = domain.manager.xmax + ymax = domain.manager.ymax + zmax = domain.manager.zmax + periodic_in_x = domain.manager.periodic_in_x + periodic_in_y = domain.manager.periodic_in_y + periodic_in_z = domain.manager.periodic_in_z + gamma_yx = domain.manager.gamma_yx + gamma_zx = domain.manager.gamma_zx + gamma_zy = domain.manager.gamma_zy + n_layers = domain.manager.n_layers + N = self.steps or int(ceil(self.dt/self.fiber_dt)) + # dt = self.dt/N + self.domain = DomainManager(xmin=xmin, xmax=xmax, ymin=ymin, + ymax=ymax, zmin=zmin, zmax=zmax, + periodic_in_x=periodic_in_x, + periodic_in_y=periodic_in_y, + periodic_in_z=periodic_in_z, + gamma_yx=gamma_yx, + gamma_zx=gamma_zx, + gamma_zy=gamma_zy, + n_layers=n_layers, + dt=self.dt, + calls_per_step=N + ) + else: + self.domain = None + # A seperate list for the nearest neighbourhood search is + # benefitial since it is much smaller than the original one. + nnps = LinkedListNNPS(dim=scheme.dim, particles=particles, + radius_scale=kernel.radius_scale, + domain=self.domain, + fixed_h=False, cache=False, sort_gids=False) + # The acceleration evaluator needs to be set up in order to compile + # it together with the integrator. + if parallel: + self.acceleration_eval = AccelerationEval( + particle_arrays=particles, + equations=equations, + kernel=kernel) + else: + self.acceleration_eval = AccelerationEval( + particle_arrays=particles, + equations=equations, + kernel=kernel, + mode='serial') + # Compilation of the integrator not using openmp, because the + # overhead is too large for those few fiber particles. + comp = SPHCompiler(self.acceleration_eval, self.fiber_integrator) + if parallel: + comp.compile() + else: + config = get_config() + config.use_openmp = False + comp.compile() + config.use_openmp = True + self.acceleration_eval.set_nnps(nnps) + + # Connecting neighbourhood list to integrator. + self.fiber_integrator.set_nnps(nnps) + + def post_stage(self, current_time, dt, stage): + """This post stage function gets called after each outer loop and + starts an inner loop for the fiber iteration.""" + from math import ceil + if self.innerloop: + # 1) predictor + # 2) post stage 1: + if stage == 1: + N = self.steps or int(ceil(self.dt/self.fiber_dt)) + for n in range(0, N): + self.fiber_integrator.step(current_time, dt/N) + current_time += dt/N + if self.domain_updates and self.domain: + self.domain.update() + # 3) Evaluation + # 4) post stage 2 + + class DensityCorrection(Tool): """ A tool to reinitialize the density of the fluid particles diff --git a/pysph/solver/utils.py b/pysph/solver/utils.py index 1de666633..de114c4be 100644 --- a/pysph/solver/utils.py +++ b/pysph/solver/utils.py @@ -1,6 +1,7 @@ """ Module contains some common functions. """ +from __future__ import unicode_literals # standard imports import errno @@ -130,7 +131,7 @@ def __init__(self, ti, tf, show=True, file=None, ascii=False): def _fmt_bar(self, percent, width): chars = ASCII_FMT if self.ascii else UTF_FMT nsyms = len(chars) - 1 - tens, ones = divmod(int(percent/100 * width * nsyms), nsyms) + tens, ones = divmod(int(percent/100.0 * width * nsyms), nsyms) end = chars[ones] if ones > 0 else '' return (chars[-1]*tens + end).ljust(width) diff --git a/pysph/solver/vtk_output.py b/pysph/solver/vtk_output.py index c47203821..11e746800 100644 --- a/pysph/solver/vtk_output.py +++ b/pysph/solver/vtk_output.py @@ -28,7 +28,7 @@ def set_output_vector(self, **vectors): vectors: Vectors to dump - Example V=['u', 'v', 'z'] + Example velocity=['u', 'v', 'z'] """ self.vectors = {} @@ -136,7 +136,7 @@ def dump_vtk(filename, particles, scalars=None, **vectors): vectors: Vectors to dump - Example V=['u', 'v', 'z'] + Example velocity=['u', 'v', 'z'] """ if has_pyvisfile(): diff --git a/pysph/sph/fiber/__init__.py b/pysph/sph/fiber/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/pysph/sph/fiber/beadchain.py b/pysph/sph/fiber/beadchain.py new file mode 100644 index 000000000..77dcab734 --- /dev/null +++ b/pysph/sph/fiber/beadchain.py @@ -0,0 +1,289 @@ +"""Equations for fibers based on bead chain model. + +Reference +--------- + + .. [Meyer2020] N. Meyer et. al "Parameter Identification of Fiber + Orientation Models Based on Direct Fiber Simulation with Smoothed + Particle Hydrodynamics", + Journal of Composites Science, 2020, 4, 77; doi:10.3390/jcs4020077 + +""" + +from math import acos, sin, sqrt +from numpy import pi as M_PI + +from pysph.sph.equation import Equation + + +class EBGVelocityReset(Equation): + """Reset EBG velocities. + + This is only necessary, if subcycling is used. + """ + + def loop(self, d_idx, d_eu, d_ev, d_ew, d_ex, d_ey, d_ez, + d_u, d_v, d_w, dt): + d_eu[d_idx] = 0 + d_ev[d_idx] = 0 + d_ew[d_idx] = 0 + + d_u[d_idx] += d_ex[d_idx]/dt + d_v[d_idx] += d_ey[d_idx]/dt + d_w[d_idx] += d_ez[d_idx]/dt + + d_ex[d_idx] = 0 + d_ey[d_idx] = 0 + d_ez[d_idx] = 0 + + +class Tension(Equation): + """Linear elastic fiber tension. + + Particle acceleration based on fiber tension is computed. The source + must be chosen to be the same as the destination particles. + See eq. (16) in [Meyer2020]. + """ + + def __init__(self, dest, sources, ea): + r""" + Parameters + ---------- + ea : float + rod stiffness (elastic modulus x section area) + """ + self.ea = ea + super(Tension, self).__init__(dest, sources) + + def initialize(self, d_idx, d_au, d_av, d_aw): + d_au[d_idx] = 0.0 + d_av[d_idx] = 0.0 + d_aw[d_idx] = 0.0 + + def loop(self, d_idx, s_idx, d_m, d_lprev, d_lnext, s_fractag, d_fractag, + d_au, d_av, d_aw, d_rxnext, d_rynext, d_rznext, d_rnext, d_rxprev, + d_ryprev, d_rzprev, d_rprev): + + # previous particle + if d_rprev[d_idx] > 1E-14: + t = self.ea*(d_rprev[d_idx]/d_lprev[d_idx]-1) + d_au[d_idx] += (t*d_rxprev[d_idx]/d_rprev[d_idx])/d_m[d_idx] + d_av[d_idx] += (t*d_ryprev[d_idx]/d_rprev[d_idx])/d_m[d_idx] + d_aw[d_idx] += (t*d_rzprev[d_idx]/d_rprev[d_idx])/d_m[d_idx] + + # next particle + if d_rnext[d_idx] > 1E-14: + t = self.ea*(d_rnext[d_idx]/d_lnext[d_idx]-1) + d_au[d_idx] += (t*d_rxnext[d_idx]/d_rnext[d_idx])/d_m[d_idx] + d_av[d_idx] += (t*d_rynext[d_idx]/d_rnext[d_idx])/d_m[d_idx] + d_aw[d_idx] += (t*d_rznext[d_idx]/d_rnext[d_idx])/d_m[d_idx] + + +class ArtificialDamping(Equation): + """Damp EBG particle motion. + + EBG Particles are damped based on EBG velocity. Use this in combination + with EBGStep and EBGVelocityReset, if subcycles are applied. + """ + + def __init__(self, dest, sources, d): + r""" + Parameters + ---------- + d : float + damping coefficient + """ + self.d = d + super(ArtificialDamping, self).__init__(dest, sources) + + def initialize(self, d_idx, d_au, d_av, d_aw): + d_au[d_idx] = 0.0 + d_av[d_idx] = 0.0 + d_aw[d_idx] = 0.0 + + def loop(self, d_idx, d_m, d_eu, d_ev, d_ew, d_au, d_av, d_aw): + d_au[d_idx] -= 2*self.d*d_eu[d_idx]/d_m[d_idx] + d_av[d_idx] -= 2*self.d*d_ev[d_idx]/d_m[d_idx] + d_aw[d_idx] -= 2*self.d*d_ew[d_idx]/d_m[d_idx] + + +class Bending(Equation): + r"""Linear elastic fiber bending + + Particle acceleration based on fiber bending is computed. The source + particles must be chosen to be the same as the destination particles. + See eq. (17) in [Meyer2020]. + """ + + def __init__(self, dest, sources, ei): + """ + Parameters + ---------- + ei : float + bending stiffness (elastic modulus x 2nd order moment) + """ + self.ei = ei + super(Bending, self).__init__(dest, sources) + + def initialize(self, d_idx, d_au, d_av, d_aw): + d_au[d_idx] = 0.0 + d_av[d_idx] = 0.0 + d_aw[d_idx] = 0.0 + + def loop(self, d_idx, d_m, d_phi0, d_fractag, d_phifrac, + d_rxnext, d_rynext, d_rznext, d_rnext, + d_rxprev, d_ryprev, d_rzprev, d_rprev, + d_au, d_av, d_aw): + if d_rnext[d_idx] > 1E-14 and d_rprev[d_idx] > 1E-14: + # vector to previous particle + xab = d_rxprev[d_idx] + yab = d_ryprev[d_idx] + zab = d_rzprev[d_idx] + rab = d_rprev[d_idx] + # vector to next particle + xbc = d_rxnext[d_idx] + ybc = d_rynext[d_idx] + zbc = d_rznext[d_idx] + rbc = d_rnext[d_idx] + + # normed dot product between vectors + # (limited to catch round off errors) + dot_prod_norm = (xab*xbc+yab*ybc+zab*zbc)/(rab*rbc) + dot_prod_norm = max(-1, dot_prod_norm) + dot_prod_norm = min(1, dot_prod_norm) + # angle between vectors + phi = acos(dot_prod_norm) + # direction of angle from cross product + norm = rab*rbc*sin(phi) + nx = (yab*zbc-zab*ybc)/norm + ny = (zab*xbc-xab*zbc)/norm + nz = (xab*ybc-yab*xbc)/norm + + # momentum + Mx = 2*self.ei*(phi-d_phi0[d_idx])/(rab+rbc)*nx + My = 2*self.ei*(phi-d_phi0[d_idx])/(rab+rbc)*ny + Mz = 2*self.ei*(phi-d_phi0[d_idx])/(rab+rbc)*nz + + if abs(phi-d_phi0[d_idx]) > d_phifrac[d_idx]: + d_fractag[d_idx] = 1 + + # forces on neighbouring particles + Fabx = (My*zab-Mz*yab)/(rab**2) + Faby = (Mz*xab-Mx*zab)/(rab**2) + Fabz = (Mx*yab-My*xab)/(rab**2) + Fbcx = (My*zbc-Mz*ybc)/(rbc**2) + Fbcy = (Mz*xbc-Mx*zbc)/(rbc**2) + Fbcz = (Mx*ybc-My*xbc)/(rbc**2) + + d_au[d_idx] -= (Fabx-Fbcx)/d_m[d_idx] + d_av[d_idx] -= (Faby-Fbcy)/d_m[d_idx] + d_aw[d_idx] -= (Fabz-Fbcz)/d_m[d_idx] + d_au[d_idx+1] -= Fbcx/d_m[d_idx+1] + d_av[d_idx+1] -= Fbcy/d_m[d_idx+1] + d_aw[d_idx+1] -= Fbcz/d_m[d_idx+1] + d_au[d_idx-1] += Fabx/d_m[d_idx-1] + d_av[d_idx-1] += Faby/d_m[d_idx-1] + d_aw[d_idx-1] += Fabz/d_m[d_idx-1] + + +class Friction(Equation): + """Fiber bending due to friction on fictive surfaces + + Since the fiber represented by a beadchain of particles has no thickness, a + term has do compensate fritction due to shear on the particles surface. + The source particles must be chosen to be the same as the destination + particles. + See Appendix A in [Meyer2020]. + """ + + def __init__(self, dest, sources, J, dx, mu, d): + r""" + Parameters + ---------- + J : float + moment of inertia + dx : float + length of segment + mu : float + absolute viscosity + d : float + fiber diameter + ar : float + fiber aspect ratio + """ + self.J = J + self.dx = dx + self.mu = mu + self.d = d + super(Friction, self).__init__(dest, sources) + + def initialize(self, d_idx, d_au, d_av, d_aw): + d_au[d_idx] = 0.0 + d_av[d_idx] = 0.0 + d_aw[d_idx] = 0.0 + + def loop(self, d_idx, d_m, d_rho, d_rxnext, d_rynext, d_rznext, + d_rnext, d_rxprev, d_ryprev, d_rzprev, d_rprev, d_fractag, + d_au, d_av, d_aw, d_dudx, d_dudy, d_dudz, d_dvdx, + d_dvdy, d_dvdz, d_dwdx, d_dwdy, d_dwdz): + if (d_rnext[d_idx] > 1E-14 and d_rprev[d_idx] > 1E-14): + + dx = d_rxprev[d_idx]-d_rxnext[d_idx] + dy = d_ryprev[d_idx]-d_rynext[d_idx] + dz = d_rzprev[d_idx]-d_rznext[d_idx] + r = sqrt(dx**2+dy**2+dz**2) + s1 = dx/r + s2 = dy/r + s3 = dz/r + + # ensuring that [sx sy sz] is not parallel to [1 0 0] + if abs(s2) > 1E-14 or abs(s3) > 1E-14: + fac = (3.*self.dx*self.d**2/4.*self.mu*M_PI)/(s2**2+s3**2) + + Mx = fac*((s1*s2**2*s3+s1*s3**3)*d_dvdx[d_idx] + + (-s1**2*s2*s3+s2*s3)*d_dvdy[d_idx] + + (-s1**2*s3**2-s2**2)*d_dvdz[d_idx] + + (-s1*s2**3-s1*s2*s3**2)*d_dwdx[d_idx] + + (s1**2*s2**2+s3**2)*d_dwdy[d_idx] + + (s1**2*s2*s3-s2*s3)*d_dwdz[d_idx]) + My = fac*((-s1*s2**2*s3-s1*s3**3)*d_dudx[d_idx] + + (s1**2*s2*s3-s2*s3)*d_dudy[d_idx] + + (s1**2*s3**2+s2**2)*d_dudz[d_idx] + + (-s2**4-2*s2**2*s3**2-s3**4)*d_dwdx[d_idx] + + (s1*s2**3+s1*s2*s3**2)*d_dwdy[d_idx] + + (s1*s2**2*s3+s1*s3**3)*d_dwdz[d_idx]) + Mz = fac*((s1*s2**3+s1*s2*s3**2)*d_dudx[d_idx] + + (-s1**2*s2**2-s3**2)*d_dudy[d_idx] + + (-s1**2*s2*s3+s2*s3)*d_dudz[d_idx] + + (s2**4+2*s2**2*s3**2+s3**4)*d_dvdx[d_idx] + + (-s1*s2**3-s1*s2*s3**2)*d_dvdy[d_idx] + + (-s1*s2**2*s3-s1*s3**3)*d_dvdz[d_idx]) + else: + fac = (3.*self.dx*self.d**2/4.*self.mu*M_PI)/(s1**2+s3**2) + + Mx = fac*((-s1*s2**2*s3+s1*s3)*d_dvdx[d_idx] + + (s1**2*s2*s3+s2*s3**3)*d_dvdy[d_idx] + + (-s2**2*s3**2-s1**2)*d_dvdz[d_idx] + + (-s1**3*s2-s1*s2*s3**2)*d_dwdx[d_idx] + + (s1**4+2*s1**2*s3**2+s3**4)*d_dwdy[d_idx] + + (-s1**2*s2*s3-s2*s3**3)*d_dwdz[d_idx]) + My = fac*((s1*s2**2*s3-s1*s3)*d_dudx[d_idx] + + (-s1**2*s2*s3-s2*s3**3)*d_dudy[d_idx] + + (s2**2*s3**2+s1**2)*d_dudz[d_idx] + + (-s1**2*s2**2-s3**2)*d_dwdx[d_idx] + + (s1**3*s2+s1*s2*s3**2)*d_dwdy[d_idx] + + (-s1*s2**2*s3+s1*s3)*d_dwdz[d_idx]) + Mz = fac*((s1**3*s2+s1*s2*s3**2)*d_dudx[d_idx] + + (-s1**4-2*s1**2*s3**2-s3**4)*d_dudy[d_idx] + + (s1**2*s2*s3+s2*s3**3)*d_dudz[d_idx] + + (s1**2*s2**2+s3**2)*d_dvdx[d_idx] + + (-s1**3*s2-s1*s2*s3**2)*d_dvdy[d_idx] + + (s1*s2**2*s3-s1*s3)*d_dvdz[d_idx]) + + d_au[d_idx+1] += (My*d_rznext[d_idx]-Mz*d_rynext[d_idx])/(2*self.J) + d_av[d_idx+1] += (Mz*d_rxnext[d_idx]-Mx*d_rznext[d_idx])/(2*self.J) + d_aw[d_idx+1] += (Mx*d_rynext[d_idx]-My*d_rxnext[d_idx])/(2*self.J) + + d_au[d_idx-1] += (My*d_rzprev[d_idx]-Mz*d_ryprev[d_idx])/(2*self.J) + d_av[d_idx-1] += (Mz*d_rxprev[d_idx]-Mx*d_rzprev[d_idx])/(2*self.J) + d_aw[d_idx-1] += (Mx*d_ryprev[d_idx]-My*d_rxprev[d_idx])/(2*self.J) diff --git a/pysph/sph/fiber/utils.py b/pysph/sph/fiber/utils.py new file mode 100644 index 000000000..94e641c98 --- /dev/null +++ b/pysph/sph/fiber/utils.py @@ -0,0 +1,541 @@ +"""Utitlity equations for fibers. + +Reference +--------- + + .. [Meyer2020] N. Meyer et. al "Parameter Identification of Fiber + Orientation Models Based on Direct Fiber Simulation with Smoothed + Particle Hydrodynamics", + Journal of Composites Science, 2020, 4, 77; doi:10.3390/jcs4020077 + +""" +from math import acos, sin, sqrt +from numpy import pi as M_PI + +from pysph.sph.equation import Equation + + +class ComputeDistance(Equation): + """Compute Distances to neighbours. + + The loop saves vectors to previous and next particle only. + """ + + def loop(self, d_idx, s_idx, d_rxnext, d_rynext, d_rznext, d_rnext, + d_rxprev, d_ryprev, d_rzprev, d_rprev, s_fractag, d_fidx, + s_fidx, d_fractag, s_rxnext, s_rynext, s_rznext, s_rnext, + s_rxprev, s_ryprev, s_rzprev, s_rprev, XIJ, RIJ): + if d_fidx[d_idx] == s_fidx[s_idx]+1: + if s_fractag[s_idx] == 0: + d_rxprev[d_idx] = -XIJ[0] + d_ryprev[d_idx] = -XIJ[1] + d_rzprev[d_idx] = -XIJ[2] + d_rprev[d_idx] = RIJ + else: + d_rxprev[d_idx] = 0.0 + d_ryprev[d_idx] = 0.0 + d_rzprev[d_idx] = 0.0 + d_rprev[d_idx] = 0.0 + s_rnext[s_idx] = 0.0 + s_rxnext[s_idx] = 0.0 + s_rynext[s_idx] = 0.0 + s_rznext[s_idx] = 0.0 + if d_fidx[d_idx] == s_fidx[s_idx]-1: + if d_fractag[d_idx] == 0: + d_rxnext[d_idx] = -XIJ[0] + d_rynext[d_idx] = -XIJ[1] + d_rznext[d_idx] = -XIJ[2] + d_rnext[d_idx] = RIJ + else: + s_rxprev[s_idx] = 0.0 + s_ryprev[s_idx] = 0.0 + s_rzprev[s_idx] = 0.0 + s_rprev[s_idx] = 0.0 + d_rxnext[d_idx] = 0.0 + d_rynext[d_idx] = 0.0 + d_rznext[d_idx] = 0.0 + d_rnext[d_idx] = 0.0 + + +class HoldPoints(Equation): + """Holds flagged points. + + Points tagged with 'holdtag' == tag are excluded from accelaration. This + little trick allows testing of fibers with fixed BCs. + """ + + def __init__(self, dest, sources, tag, x=True, y=True, z=True, + mirror_particle=0): + r""" + Parameters + ---------- + tags : int + tag of fixed particle defined as property 'holdtag' + x : boolean + True, if x-position should not be changed + y : boolean + True, if y-position should not be changed + z : boolean + True, if z-position should not be changed + mirror_particle : int + idx shift to a particle, which displacement should be mirrored to + origin + """ + self.tag = tag + self.x = x + self.y = y + self.z = z + self.mirror = mirror_particle + super(HoldPoints, self).__init__(dest, sources) + + def loop(self, d_idx, d_holdtag, d_au, d_av, d_aw, d_auhat, d_avhat, + d_awhat, d_u, d_v, d_w, d_x, d_y, d_z, d_Fx, d_Fy, d_Fz, d_m): + if d_holdtag[d_idx] == self.tag: + if self.x: + d_Fx[d_idx] = d_m[d_idx] * d_au[d_idx] + d_au[d_idx] = 0 + d_auhat[d_idx] = 0 + d_u[d_idx] = 0 + elif not self.mirror == 0: + # Copy properties to mirror particle + d_x[d_idx+self.mirror] = -d_x[d_idx] + d_u[d_idx+self.mirror] = -d_u[d_idx] + + if self.y: + d_Fy[d_idx] = d_m[d_idx] * d_av[d_idx] + d_av[d_idx] = 0 + d_avhat[d_idx] = 0 + d_v[d_idx] = 0 + elif not self.mirror == 0: + # Copy properties to mirror particle + d_y[d_idx+self.mirror] = -d_y[d_idx] + d_v[d_idx+self.mirror] = -d_v[d_idx] + + if self.z: + d_Fz[d_idx] = d_m[d_idx] * d_aw[d_idx] + d_aw[d_idx] = 0 + d_awhat[d_idx] = 0 + d_w[d_idx] = 0 + elif not self.mirror == 0: + # Copy properties to mirror particle + d_z[d_idx+self.mirror] = -d_z[d_idx] + d_w[d_idx+self.mirror] = -d_w[d_idx] + + +class Vorticity(Equation): + """Computes vorticity of velocity field + + According to Monaghan 1992 (2.12). + """ + + def initialize(self, d_idx, d_omegax, d_omegay, d_omegaz): + d_omegax[d_idx] = 0.0 + d_omegay[d_idx] = 0.0 + d_omegaz[d_idx] = 0.0 + + def loop(self, d_idx, s_idx, d_rho, s_m, d_omegax, d_omegay, d_omegaz, + DWIJ, VIJ): + v = s_m[s_idx]/d_rho[d_idx] + d_omegax[d_idx] += v*(VIJ[1]*DWIJ[2]-VIJ[2]*DWIJ[1]) + d_omegay[d_idx] += v*(VIJ[2]*DWIJ[0]-VIJ[0]*DWIJ[2]) + d_omegaz[d_idx] += v*(VIJ[0]*DWIJ[1]-VIJ[1]*DWIJ[0]) + + +class VelocityGradient(Equation): + """Computes 2nd order tensor representing the velocity gradient. + + See eq. (25) in [Meyer2020]. + """ + + def initialize(self, d_idx, d_dudx, d_dudy, d_dudz, d_dvdx, d_dvdy, d_dvdz, + d_dwdx, d_dwdy, d_dwdz): + d_dudx[d_idx] = 0.0 + d_dudy[d_idx] = 0.0 + d_dudz[d_idx] = 0.0 + + d_dvdx[d_idx] = 0.0 + d_dvdy[d_idx] = 0.0 + d_dvdz[d_idx] = 0.0 + + d_dwdx[d_idx] = 0.0 + d_dwdy[d_idx] = 0.0 + d_dwdz[d_idx] = 0.0 + + def loop(self, d_idx, s_idx, s_rho, s_m, d_dudx, d_dudy, d_dudz, d_dvdx, + d_dvdy, d_dvdz, d_dwdx, d_dwdy, d_dwdz, DWIJ, VIJ): + v = s_m[s_idx]/s_rho[s_idx] + d_dudx[d_idx] -= v*VIJ[0]*DWIJ[0] + d_dudy[d_idx] -= v*VIJ[0]*DWIJ[1] + d_dudz[d_idx] -= v*VIJ[0]*DWIJ[2] + + d_dvdx[d_idx] -= v*VIJ[1]*DWIJ[0] + d_dvdy[d_idx] -= v*VIJ[1]*DWIJ[1] + d_dvdz[d_idx] -= v*VIJ[1]*DWIJ[2] + + d_dwdx[d_idx] -= v*VIJ[2]*DWIJ[0] + d_dwdy[d_idx] -= v*VIJ[2]*DWIJ[1] + d_dwdz[d_idx] -= v*VIJ[2]*DWIJ[2] + + +class Damping(Equation): + """Damp particle motion. + + Particles are damped. Difference to ArtificialDamping: This damps real + particle velocities and therefore affects not only the fiber iteration. + In this context it may be used to test fiber contact in a damped + environment. + """ + + def __init__(self, dest, sources, d): + r""" + Parameters + ---------- + d : float + damping coefficient + """ + self.d = d + super(Damping, self).__init__(dest, sources) + + def initialize(self, d_idx, d_au, d_av, d_aw): + d_au[d_idx] = 0.0 + d_av[d_idx] = 0.0 + d_aw[d_idx] = 0.0 + + def loop(self, d_idx, d_m, d_u, d_v, d_w, d_au, d_av, d_aw): + d_au[d_idx] -= 2*self.d*d_u[d_idx]/d_m[d_idx] + d_av[d_idx] -= 2*self.d*d_v[d_idx]/d_m[d_idx] + d_aw[d_idx] -= 2*self.d*d_w[d_idx]/d_m[d_idx] + + +class SimpleContact(Equation): + """This class computes simple fiber repulsion to stop penetration. + + It computes the force between two spheres as Hertz pressure. + """ + + def __init__(self, dest, sources, E, d, pois=0.3): + r""" + Parameters + ---------- + E : float + Young's modulus + d : float + fiber diameter + pois : flost + poisson number + """ + self.E = E + self.d = d + self.pois = pois + super(SimpleContact, self).__init__(dest, sources) + + def initialize(self, d_idx, d_au, d_av, d_aw): + d_au[d_idx] = 0.0 + d_av[d_idx] = 0.0 + d_aw[d_idx] = 0.0 + + def loop(self, d_idx, s_idx, d_fidx, s_fidx, d_m, d_au, d_av, d_aw, + XIJ, RIJ): + if not s_fidx[s_idx] == d_fidx[d_idx] and RIJ < self.d: + E_star = 1/(2*((1-self.pois**2)/self.E)) + # effective radius for two spheres of same size + R = self.d/4 + F = 4/3 * E_star * sqrt(R) * abs(self.d-RIJ)**1.5 + d_au[d_idx] += XIJ[0]/RIJ * F/d_m[d_idx] + d_av[d_idx] += XIJ[1]/RIJ * F/d_m[d_idx] + d_aw[d_idx] += XIJ[2]/RIJ * F/d_m[d_idx] + + +class Contact(Equation): + """This class computes fiber repulsion to stop penetration. + + It computes the force between two spheres based on Hertz pressure + between two cylinders. This Equation requires a computation of + distances by the Bending equation. + + See eq. (27)-(34) in [Meyer2020]. + """ + def __init__(self, dest, sources, E, d, dim, pois=0.3, k=0.0, lim=0.1, + eta0=0.0, dt=0.0): + r""" + Parameters + ---------- + E : float + Young's modulus + d : float + fiber diameter + dim : int + dimensionaltiy of the problem + pois : float + poisson number + k : float + friction coefficient between fibers + eta0 : float + viscosity of suspension fluid + """ + self.E = E + self.d = d + self.pois = pois + self.k = k + self.dim = dim + self.lim = lim + self.eta0 = eta0 + self.dt = dt + self.Fx = 0.0 + self.Fy = 0.0 + self.Fz = 0.0 + self.nx = 0.0 + self.ny = 0.0 + self.nz = 0.0 + self.proj = 0.0 + self.dist = 0.0 + super(Contact, self).__init__(dest, sources) + + def initialize(self, d_idx, d_au, d_av, d_aw, d_Fx, d_Fy, d_Fz): + d_au[d_idx] = 0.0 + d_av[d_idx] = 0.0 + d_aw[d_idx] = 0.0 + d_Fx[d_idx] = 0.0 + d_Fy[d_idx] = 0.0 + d_Fz[d_idx] = 0.0 + + def loop(self, d_idx, s_idx, d_m, d_au, d_av, d_aw, d_rxnext, d_rynext, + d_rznext, d_rnext, d_rxprev, d_ryprev, d_rzprev, d_rprev, + s_rxnext, s_rynext, s_rznext, s_rnext, s_rxprev, s_ryprev, + s_rzprev, s_rprev, d_Fx, d_Fy, d_Fz, d_fractag, d_tag, s_tag, + XIJ, VIJ, RIJ, EPS): + """Loop over all particles and compute interaction forces.""" + + # exclude self and direct neighbors + if (RIJ > 1E-14 + and RIJ < 1.5*self.d + and abs(RIJ-d_rprev[d_idx]) > 1E-14 + and abs(RIJ-d_rnext[d_idx]) > 1E-14): + + # unit direction of destination fiber + dx = d_rxprev[d_idx] - d_rxnext[d_idx] + dy = d_ryprev[d_idx] - d_rynext[d_idx] + dz = d_rzprev[d_idx] - d_rznext[d_idx] + dr = sqrt(dx**2 + dy**2 + dz**2) + dx = dx/dr + dy = dy/dr + dz = dz/dr + + # unit direction of source fiber + sx = s_rxprev[s_idx] - s_rxnext[s_idx] + sy = s_ryprev[s_idx] - s_rynext[s_idx] + sz = s_rzprev[s_idx] - s_rznext[s_idx] + sr = sqrt(sx**2 + sy**2 + sz**2) + sx = sx/sr + sy = sy/sr + sz = sz/sr + + # conditions + d_prev_tip = (d_rprev[d_idx] < EPS + and dx*XIJ[0]+dy*XIJ[1]+dz*XIJ[2] < EPS) + s_prev_tip = (s_rprev[s_idx] < EPS + and sx*XIJ[0]+sy*XIJ[1]+sz*XIJ[2] > EPS) + d_next_tip = (d_rnext[d_idx] < EPS + and dx*XIJ[0]+dy*XIJ[1]+dz*XIJ[2] > EPS) + s_next_tip = (s_rnext[s_idx] < EPS + and sx*XIJ[0]+sy*XIJ[1]+sz*XIJ[2] < EPS) + + # determine case + if d_prev_tip or d_next_tip or dr < EPS: + if s_prev_tip or s_next_tip or sr < EPS: + alpha = self.compute_point_point_props(XIJ, RIJ) + w = 1.0 + else: + alpha = self.compute_center_point_props(XIJ, sx, sy, sz) + alpha = 0.0 + w = self.weight(s_rnext[s_idx], s_rprev[s_idx], self.proj) + elif s_prev_tip or s_next_tip or sr < EPS: + if d_prev_tip or d_next_tip or dr < EPS: + alpha = self.compute_point_point_props(XIJ, RIJ) + w = 1.0 + else: + alpha = self.compute_center_point_props(XIJ, dx, dy, dz) + w = self.weight(d_rnext[d_idx], d_rprev[d_idx], self.proj) + else: + # center and center + alpha = self.compute_center_center_props( + XIJ, dx, dy, dz, sx, sy, sz) + w = self.weight(d_rnext[d_idx], d_rprev[d_idx], self.proj) + + d = self.d - self.dist + if d >= 0.0: + self.compute_contact_force(VIJ[0], VIJ[1], VIJ[2], w, d) + elif alpha > 0.0: + pass + # self.compute_lubrication_force( + # VIJ[0], VIJ[1], VIJ[2], w, d, alpha) + else: + self.Fx = 0.0 + self.Fy = 0.0 + self.Fz = 0.0 + + d_Fx[d_idx] += self.Fx + d_Fy[d_idx] += self.Fy + d_Fz[d_idx] += self.Fz + + d_au[d_idx] += self.Fx/d_m[d_idx] + d_av[d_idx] += self.Fy/d_m[d_idx] + d_aw[d_idx] += self.Fz/d_m[d_idx] + + def compute_center_point_props(self, XIJ, ax=0.0, ay=0.0, az=0.0): + """Determine normal, distance and projected length - Center-Point. + + The projection looks like this: + + x (source point ) + /| + XIJ/ | distance + / | + dest point o---*----> unit vector a + proj + """ + projection = -(ax*XIJ[0] + ay*XIJ[1] + az*XIJ[2]) + tx = XIJ[0] + projection*ax + ty = XIJ[1] + projection*ay + tz = XIJ[2] + projection*az + tr = sqrt(tx**2 + ty**2 + tz**2) + + self.nx = tx/tr + self.ny = ty/tr + self.nz = tz/tr + self.proj = projection + self.dist = tr + + return 0.0 + + def compute_point_point_props(self, XIJ, RIJ): + """Compute the normal between end points. + + This case is trivial. + """ + self.nx = XIJ[0]/RIJ + self.ny = XIJ[1]/RIJ + self.nz = XIJ[2]/RIJ + self.dist = RIJ + + return 0.0 + + def compute_center_center_props(self, XIJ, dx=0.0, dy=0.0, dz=0.0, + sx=0.0, sy=0.0, sz=0.0): + """Compute normal direction between two lines. + + The normal is given either by the cross product or by a projection, if + both lines are (almost) parallel. + """ + # normal direction at contact + nx = dy*sz - dz*sy + ny = dz*sx - dx*sz + nz = dx*sy - dy*sx + nr = sqrt(nx**2 + ny**2 + nz**2) + + # Parallel vectors should be treated with a projection + if nr > 0.01: + self.nx = nx/nr + self.ny = ny/nr + self.nz = nz/nr + alpha = acos(dx*sx+dy*sy+dz*sz) + self.project_center_center(XIJ, dx, dy, dz, sx, sy, sz) + else: + self.compute_center_point_props(XIJ, dx, dy, dz) + alpha = M_PI/2. + + return alpha + + def project_center_center(self, XIJ, dx=0.0, dy=0.0, dz=0.0, + sx=0.0, sy=0.0, sz=0.0): + """Find distance and projection of contact point on fiber line. + + Therefore, solve + |dx nx -sx| |dest projection | |XIJ[0]| + |dy ny -sy| |distance | = |XIJ[1]| + |dz nz -sz| |src projection | |XIJ[2]| + + """ + nx = self.nx + ny = self.ny + nz = self.nz + det = (dx*ny*(-sz) + nx*(-sy)*dz + (-sx)*dy*nz + - (-sx)*ny*dz - dx*(-sy)*nz - nx*dy*(-sz)) + self.proj = ((sy*nz - ny*sz)*XIJ[0] + + (nx*sz - sx*nz)*XIJ[1] + + (sx*ny - nx*sy)*XIJ[2])/det + self.dist = ((dy*sz - sy*dz)*XIJ[0] + + (sx*dz - dx*sz)*XIJ[1] + + (dx*sy - sx*dy)*XIJ[2])/det + if self.dist < 0.0: + self.dist = -self.dist + self.nx = -nx + self.ny = -ny + self.nz = -nz + + def weight(self, next=0.0, prev=0.0, proj=0.0): + """Weight force contribution among neighboring particles. + + Any projection outside the range between the previous and next particle + is weighted with zero. Anything between is linearly distributed with + special treating for the edge cases of fiber tips. + """ + w = 0.0 + if next < 1E-14 and proj > -prev and proj <= 0: + w = (prev+proj)/prev + elif prev < 1E-14 and proj < next and proj >= 0: + w = (next-proj)/next + elif proj < prev and proj >= 0: + w = (prev-proj)/prev + elif proj > -next and proj <= 0: + w = (next+proj)/next + return w + + def compute_contact_force(self, vx=0.0, vy=0.0, vz=0.0, w=0.0, d=1.0): + """Compute the contact force at interaction point.""" + + v_dot_n = vx*self.nx + vy*self.ny + vz*self.nz + + # elastic factor from Hertz' pressure in contact + E_star = 1/(2*((1-self.pois**2)/self.E)) + + d = min(self.lim*self.d, d) + + if self.dim == 3: + F = 4/3 * d**1.5 * sqrt(self.d/2) * E_star + else: + # workaround for 2D contact (2 reactangular surfaces) + F = self.E*d + + vrx = vx - v_dot_n*self.nx + vry = vy - v_dot_n*self.ny + vrz = vz - v_dot_n*self.nz + v_rel = sqrt(vrx**2 + vry**2 + vrz**2) + v_rel = v_rel if v_rel > 1E-14 else 1 + + self.Fx = w*(F*self.nx - self.k*F*vrx/v_rel) + self.Fy = w*(F*self.ny - self.k*F*vry/v_rel) + self.Fz = w*(F*self.nz - self.k*F*vrz/v_rel) + + def compute_lubrication_force(self, vx=0.0, vy=0.0, vz=0.0, w=0.0, d=1.0, + alpha=M_PI/2): + """Compute the lubrication force at interaction point.""" + # limit extreme forces + d = min(d, -0.01*self.d) + R = self.d/2 + + v_dot_n = vx*self.nx + vy*self.ny + vz*self.nz + + if abs(alpha) > 0.1: + A = R**2/abs(sin(alpha)) + F = 12.*A*M_PI*self.eta0*v_dot_n/d + else: + # Lindstroem limit to treat parallel cylinders (singularity!) + A0 = 3.*M_PI*sqrt(2.)/8. + A1 = 207*M_PI*sqrt(2.)/160. + L = self.d + F = -L*self.eta0*v_dot_n*(A0-A1*d/R)*(-d/R)**(-3./2.) + + self.Fx = w*F*self.nx + self.Fy = w*F*self.ny + self.Fz = w*F*self.nz diff --git a/pysph/sph/integrator.py b/pysph/sph/integrator.py index 1656317b6..bb1666325 100644 --- a/pysph/sph/integrator.py +++ b/pysph/sph/integrator.py @@ -6,8 +6,8 @@ from the `sph_eval` module. """ -from numpy import sqrt import numpy as np +from numpy import sqrt # Local imports. from .integrator_step import IntegratorStep diff --git a/pysph/sph/integrator_step.py b/pysph/sph/integrator_step.py index b757d4512..81abed917 100644 --- a/pysph/sph/integrator_step.py +++ b/pysph/sph/integrator_step.py @@ -11,8 +11,9 @@ class IntegratorStep(object): """Subclass this and implement the methods ``initialize``, ``stage1`` etc. Use the same conventions as the equations. """ + def __repr__(self): - return '%s()'%(self.__class__.__name__) + return '%s()' % (self.__class__.__name__) ############################################################################### @@ -21,7 +22,7 @@ def __repr__(self): class EulerStep(IntegratorStep): """Fast but inaccurate integrator. Use this for testing""" def stage1(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_x, d_y, - d_z, d_rho, d_arho, dt): + d_z, d_rho, d_arho, dt): d_u[d_idx] += dt*d_au[d_idx] d_v[d_idx] += dt*d_av[d_idx] d_w[d_idx] += dt*d_aw[d_idx] @@ -32,6 +33,77 @@ def stage1(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_x, d_y, d_rho[d_idx] += dt*d_arho[d_idx] + +############################################################################### +# `EBGStep` class +############################################################################### +class EBGStep(IntegratorStep): + """Euler integrator using alternative velocity for EBG.""" + def stage1(self, d_idx, d_eu, d_ev, d_ew, d_au, d_av, d_aw, d_x, d_y, + d_z, d_ex, d_ey, d_ez, dt): + d_eu[d_idx] += dt*d_au[d_idx] + d_ev[d_idx] += dt*d_av[d_idx] + d_ew[d_idx] += dt*d_aw[d_idx] + + d_x[d_idx] += dt*d_eu[d_idx] + d_y[d_idx] += dt*d_ev[d_idx] + d_z[d_idx] += dt*d_ew[d_idx] + + d_ex[d_idx] += dt*d_eu[d_idx] + d_ey[d_idx] += dt*d_ev[d_idx] + d_ez[d_idx] += dt*d_ew[d_idx] + + +############################################################################### +# `RK2Step` class +############################################################################### +class RK2Step(IntegratorStep): + r"""General 2nd order Runge-Kutta Integrator + + In the predictor step the particles are advanced to `t + dt/2`. The + particles are then advanced with the new force computed at this position. + .. math:: + + y_{n+\frac{1}{2}} = y_n + \frac{\Delta t}{2} f(t_n, y_n) + y_{n+1} = y_n + \Delta t f{y_{n+\frac{1}{2}}} + + This integrator should be used in EPEC mode to mimic a 2nd order + Runge-Kutta Scheme. + + """ + def initialize(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, + d_u0, d_v0, d_w0, d_u, d_v, d_w): + d_x0[d_idx] = d_x[d_idx] + d_y0[d_idx] = d_y[d_idx] + d_z0[d_idx] = d_z[d_idx] + + d_u0[d_idx] = d_u[d_idx] + d_v0[d_idx] = d_v[d_idx] + d_w0[d_idx] = d_w[d_idx] + + def stage1(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, + d_u0, d_v0, d_w0, d_u, d_v, d_w, d_au, d_av, d_aw, dt): + dtb2 = 0.5*dt + d_u[d_idx] = d_u0[d_idx] + dtb2*d_au[d_idx] + d_v[d_idx] = d_v0[d_idx] + dtb2*d_av[d_idx] + d_w[d_idx] = d_w0[d_idx] + dtb2*d_aw[d_idx] + + d_x[d_idx] = d_x0[d_idx] + dtb2 * d_u[d_idx] + d_y[d_idx] = d_y0[d_idx] + dtb2 * d_v[d_idx] + d_z[d_idx] = d_z0[d_idx] + dtb2 * d_w[d_idx] + + def stage2(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, + d_u0, d_v0, d_w0, d_u, d_v, d_w, d_au, d_av, d_aw, dt): + + d_u[d_idx] = d_u0[d_idx] + dt * d_au[d_idx] + d_v[d_idx] = d_v0[d_idx] + dt * d_av[d_idx] + d_w[d_idx] = d_w0[d_idx] + dt * d_aw[d_idx] + + d_x[d_idx] = d_x0[d_idx] + dt * d_u[d_idx] + d_y[d_idx] = d_y0[d_idx] + dt * d_v[d_idx] + d_z[d_idx] = d_z0[d_idx] + dt * d_w[d_idx] + + ############################################################################### # `WCSPHStep` class ############################################################################### @@ -61,8 +133,8 @@ def initialize(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, d_rho0[d_idx] = d_rho[d_idx] def stage1(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, - d_u0, d_v0, d_w0, d_u, d_v, d_w, d_rho0, d_rho, d_au, d_av, - d_aw, d_ax, d_ay, d_az, d_arho, dt): + d_u0, d_v0, d_w0, d_u, d_v, d_w, d_rho0, d_rho, d_au, d_av, + d_aw, d_ax, d_ay, d_az, d_arho, dt): dtb2 = 0.5*dt d_u[d_idx] = d_u0[d_idx] + dtb2*d_au[d_idx] d_v[d_idx] = d_v0[d_idx] + dtb2*d_av[d_idx] @@ -76,8 +148,8 @@ def stage1(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, d_rho[d_idx] = d_rho0[d_idx] + dtb2 * d_arho[d_idx] def stage2(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, - d_u0, d_v0, d_w0, d_u, d_v, d_w, d_rho0, d_rho, d_au, d_av, - d_aw, d_ax, d_ay, d_az, d_arho, dt): + d_u0, d_v0, d_w0, d_u, d_v, d_w, d_rho0, d_rho, d_au, d_av, + d_aw, d_ax, d_ay, d_az, d_arho, dt): d_u[d_idx] = d_u0[d_idx] + dt*d_au[d_idx] d_v[d_idx] = d_v0[d_idx] + dt*d_av[d_idx] @@ -90,6 +162,7 @@ def stage2(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, # Update densities and smoothing lengths from the accelerations d_rho[d_idx] = d_rho0[d_idx] + dt * d_arho[d_idx] + ############################################################################### # `WCSPHTVDRK3` Integrator ############################################################################### @@ -135,17 +208,18 @@ def stage2(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, d_aw, d_ax, d_ay, d_az, d_arho, dt): # update velocities - d_u[d_idx] = 0.75*d_u0[d_idx] + 0.25*( d_u[d_idx] + dt * d_au[d_idx] ) - d_v[d_idx] = 0.75*d_v0[d_idx] + 0.25*( d_v[d_idx] + dt * d_av[d_idx] ) - d_w[d_idx] = 0.75*d_w0[d_idx] + 0.25*( d_w[d_idx] + dt * d_aw[d_idx] ) + d_u[d_idx] = 0.75*d_u0[d_idx] + 0.25*(d_u[d_idx] + dt * d_au[d_idx]) + d_v[d_idx] = 0.75*d_v0[d_idx] + 0.25*(d_v[d_idx] + dt * d_av[d_idx]) + d_w[d_idx] = 0.75*d_w0[d_idx] + 0.25*(d_w[d_idx] + dt * d_aw[d_idx]) # update positions - d_x[d_idx] = 0.75*d_x0[d_idx] + 0.25*( d_x[d_idx] + dt * d_ax[d_idx] ) - d_y[d_idx] = 0.75*d_y0[d_idx] + 0.25*( d_y[d_idx] + dt * d_ay[d_idx] ) - d_z[d_idx] = 0.75*d_z0[d_idx] + 0.25*( d_z[d_idx] + dt * d_az[d_idx] ) + d_x[d_idx] = 0.75*d_x0[d_idx] + 0.25*(d_x[d_idx] + dt * d_ax[d_idx]) + d_y[d_idx] = 0.75*d_y0[d_idx] + 0.25*(d_y[d_idx] + dt * d_ay[d_idx]) + d_z[d_idx] = 0.75*d_z0[d_idx] + 0.25*(d_z[d_idx] + dt * d_az[d_idx]) # Update density - d_rho[d_idx] = 0.75*d_rho0[d_idx] + 0.25*( d_rho[d_idx] + dt * d_arho[d_idx] ) + d_rho[d_idx] = 0.75*d_rho0[d_idx] + 0.25*( + d_rho[d_idx] + dt * d_arho[d_idx]) def stage3(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, d_u0, d_v0, d_w0, d_u, d_v, d_w, d_rho0, d_rho, d_au, d_av, @@ -155,17 +229,19 @@ def stage3(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, twoby3 = 2./3. # update velocities - d_u[d_idx] = oneby3*d_u0[d_idx] + twoby3*( d_u[d_idx] + dt * d_au[d_idx] ) - d_v[d_idx] = oneby3*d_v0[d_idx] + twoby3*( d_v[d_idx] + dt * d_av[d_idx] ) - d_w[d_idx] = oneby3*d_w0[d_idx] + twoby3*( d_w[d_idx] + dt * d_aw[d_idx] ) + d_u[d_idx] = oneby3*d_u0[d_idx] + twoby3*(d_u[d_idx] + dt*d_au[d_idx]) + d_v[d_idx] = oneby3*d_v0[d_idx] + twoby3*(d_v[d_idx] + dt*d_av[d_idx]) + d_w[d_idx] = oneby3*d_w0[d_idx] + twoby3*(d_w[d_idx] + dt*d_aw[d_idx]) # update positions - d_x[d_idx] = oneby3*d_x0[d_idx] + twoby3*( d_x[d_idx] + dt * d_ax[d_idx] ) - d_y[d_idx] = oneby3*d_y0[d_idx] + twoby3*( d_y[d_idx] + dt * d_ay[d_idx] ) - d_z[d_idx] = oneby3*d_z0[d_idx] + twoby3*( d_z[d_idx] + dt * d_az[d_idx] ) + d_x[d_idx] = oneby3*d_x0[d_idx] + twoby3*(d_x[d_idx] + dt*d_ax[d_idx]) + d_y[d_idx] = oneby3*d_y0[d_idx] + twoby3*(d_y[d_idx] + dt*d_ay[d_idx]) + d_z[d_idx] = oneby3*d_z0[d_idx] + twoby3*(d_z[d_idx] + dt*d_az[d_idx]) # Update density - d_rho[d_idx] = oneby3*d_rho0[d_idx] + twoby3*( d_rho[d_idx] + dt * d_arho[d_idx] ) + d_rho[d_idx] = oneby3*d_rho0[d_idx] + twoby3*( + d_rho[d_idx] + dt*d_arho[d_idx]) + ############################################################################### # `SolidMechStep` class @@ -196,12 +272,12 @@ def initialize(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, d_s220[d_idx] = d_s22[d_idx] def stage1(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, - d_u0, d_v0, d_w0, d_u, d_v, d_w, d_rho0, d_rho, d_au, d_av, - d_aw, d_ax, d_ay, d_az, d_arho, d_e, d_e0, d_ae, - d_s00, d_s01, d_s02, d_s11, d_s12, d_s22, - d_s000, d_s010, d_s020, d_s110, d_s120, d_s220, - d_as00, d_as01, d_as02, d_as11, d_as12, d_as22, - dt): + d_u0, d_v0, d_w0, d_u, d_v, d_w, d_rho0, d_rho, d_au, d_av, + d_aw, d_ax, d_ay, d_az, d_arho, d_e, d_e0, d_ae, + d_s00, d_s01, d_s02, d_s11, d_s12, d_s22, + d_s000, d_s010, d_s020, d_s110, d_s120, d_s220, + d_as00, d_as01, d_as02, d_as11, d_as12, d_as22, + dt): dtb2 = 0.5*dt d_u[d_idx] = d_u0[d_idx] + dtb2*d_au[d_idx] d_v[d_idx] = d_v0[d_idx] + dtb2*d_av[d_idx] @@ -224,12 +300,12 @@ def stage1(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, d_s22[d_idx] = d_s220[d_idx] + dtb2 * d_as22[d_idx] def stage2(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, - d_u0, d_v0, d_w0, d_u, d_v, d_w, d_rho0, d_rho, d_au, d_av, - d_aw, d_ax, d_ay, d_az, d_arho, d_e, d_ae, d_e0, - d_s00, d_s01, d_s02, d_s11, d_s12, d_s22, - d_s000, d_s010, d_s020, d_s110, d_s120, d_s220, - d_as00, d_as01, d_as02, d_as11, d_as12, d_as22, - dt): + d_u0, d_v0, d_w0, d_u, d_v, d_w, d_rho0, d_rho, d_au, d_av, + d_aw, d_ax, d_ay, d_az, d_arho, d_e, d_ae, d_e0, + d_s00, d_s01, d_s02, d_s11, d_s12, d_s22, + d_s000, d_s010, d_s020, d_s110, d_s120, d_s220, + d_as00, d_as01, d_as02, d_as11, d_as12, d_as22, + dt): d_u[d_idx] = d_u0[d_idx] + dt*d_au[d_idx] d_v[d_idx] = d_v0[d_idx] + dt*d_av[d_idx] @@ -251,6 +327,7 @@ def stage2(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, d_s12[d_idx] = d_s120[d_idx] + dt * d_as12[d_idx] d_s22[d_idx] = d_s220[d_idx] + dt * d_as22[d_idx] + ############################################################################### # `TransportVelocityStep` class ############################################################################### @@ -265,8 +342,8 @@ class TransportVelocityStep(IntegratorStep): def initialize(self): pass - def stage1(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_uhat, d_auhat, d_vhat, - d_avhat, d_what, d_awhat, d_x, d_y, d_z, dt): + def stage1(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_uhat, d_auhat, + d_vhat, d_avhat, d_what, d_awhat, d_x, d_y, d_z, dt): dtb2 = 0.5*dt # velocity update eqn (14) @@ -296,6 +373,7 @@ def stage2(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_vmag2, dt): d_vmag2[d_idx] = (d_u[d_idx]*d_u[d_idx] + d_v[d_idx]*d_v[d_idx] + d_w[d_idx]*d_w[d_idx]) + ############################################################################### # `AdamiVerletStep` class ############################################################################### @@ -311,7 +389,8 @@ class AdamiVerletStep(IntegratorStep): def initialize(self): pass - def stage1(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_x, d_y, d_z, dt): + def stage1(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, + d_x, d_y, d_z, dt): dtb2 = 0.5*dt # velocity predictor eqn (14) @@ -345,6 +424,7 @@ def stage2(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_x, d_y, d_z, d_vmag2[d_idx] = (d_u[d_idx]*d_u[d_idx] + d_v[d_idx]*d_v[d_idx] + d_w[d_idx]*d_w[d_idx]) + ############################################################################### # `GasDFluidStep` class ############################################################################### @@ -466,7 +546,6 @@ def initialize(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, d_rho0[d_idx] = d_rho[d_idx] - def stage1(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, d_u0, d_v0, d_w0, d_u, d_v, d_w, d_e0, d_e, d_au, d_av, d_aw, d_ae, d_rho, d_rho0, d_arho, dt): @@ -499,7 +578,6 @@ def stage2(self, d_idx, d_x0, d_y0, d_z0, d_x, d_y, d_z, d_e[d_idx] = d_e0[d_idx] + dt * d_ae[d_idx] - ############################################################################### # `TwoStageRigidBodyStep` class ############################################################################### @@ -553,6 +631,7 @@ def stage2(self, d_idx, d_x, d_y, d_z, d_x0, d_y0, d_z0, d_y[d_idx] = d_y0[d_idx] + dt * 0.5 * (d_v[d_idx] + d_v0[d_idx]) d_z[d_idx] = d_z0[d_idx] + dt * 0.5 * (d_w[d_idx] + d_w0[d_idx]) + ############################################################################### # `OneStageRigidBodyStep` class ############################################################################### @@ -570,13 +649,13 @@ def initialize(self, d_idx, d_x, d_y, d_z, d_x0, d_y0, d_z0, d_z0[d_idx] = d_z[d_idx] def stage1(self, d_idx, d_x, d_y, d_z, d_x0, d_y0, d_z0, - d_u, d_v, d_w, d_u0, d_v0, d_w0, d_au, d_av, d_aw, - dt): + d_u, d_v, d_w, d_u0, d_v0, d_w0, d_au, d_av, d_aw, + dt): pass def stage2(self, d_idx, d_x, d_y, d_z, d_x0, d_y0, d_z0, - d_u, d_v, d_w, d_u0, d_v0, d_w0, d_au, d_av, d_aw, - dt): + d_u, d_v, d_w, d_u0, d_v0, d_w0, d_au, d_av, d_aw, + dt): # update velocities d_u[d_idx] += dt * d_au[d_idx] @@ -681,6 +760,7 @@ def stage2(self, d_idx, d_x, d_y, d_z, d_u, d_v, d_w, d_v[d_idx] += dtb2 * d_av[d_idx] d_w[d_idx] += dtb2 * d_aw[d_idx] + ############################################################################### # `InletOutletStep` class ############################################################################### @@ -703,14 +783,13 @@ def stage2(self, d_idx, d_x, d_y, d_z, d_u, d_v, d_w, dt): d_z[d_idx] += dtb2 * d_w[d_idx] - ############################################################################### class LeapFrogStep(IntegratorStep): r"""Using this stepper with XSPH as implemented in `pysph.base.basic_equations.XSPHCorrection` is not directly possible and - requires a nicer implementation where the correction alone is added to ``ax, - ay, az``. + requires a nicer implementation where the correction alone is added to + ``ax, ay, az``. """ def stage1(self, d_idx, d_x, d_y, d_z, d_u, d_v, d_w, d_ax, d_ay, d_az, @@ -739,8 +818,8 @@ class PEFRLStep(IntegratorStep): r"""Using this stepper with XSPH as implemented in `pysph.base.basic_equations.XSPHCorrection` is not directly possible and - requires a nicer implementation where the correction alone is added to ``ax, - ay, az``. + requires a nicer implementation where the correction alone is added to + ``ax, ay, az``. """ def stage1(self, d_idx, d_x, d_y, d_z, d_u, d_v, d_w, d_ax, d_ay, diff --git a/pysph/sph/rigid_body.py b/pysph/sph/rigid_body.py index d363998fe..bdec8db0a 100644 --- a/pysph/sph/rigid_body.py +++ b/pysph/sph/rigid_body.py @@ -303,9 +303,9 @@ def loop(self, d_idx, d_m, d_au, d_av, d_aw, d_rho, d_au[d_idx] += ax d_av[d_idx] += ay d_aw[d_idx] += az - s_fx[s_idx] += -d_m[d_idx]*ax - s_fy[s_idx] += -d_m[d_idx]*ay - s_fz[s_idx] += -d_m[d_idx]*az + s_fx[s_idx] -= d_m[d_idx]*ax + s_fy[s_idx] -= d_m[d_idx]*ay + s_fz[s_idx] -= d_m[d_idx]*az class PressureRigidBody(Equation): @@ -330,9 +330,9 @@ def loop(self, d_idx, d_m, d_rho, d_au, d_av, d_aw, d_p, d_au[d_idx] += ax d_av[d_idx] += ay d_aw[d_idx] += az - s_fx[s_idx] += -d_m[d_idx]*ax - s_fy[s_idx] += -d_m[d_idx]*ay - s_fz[s_idx] += -d_m[d_idx]*az + s_fx[s_idx] -= d_m[d_idx]*ax + s_fy[s_idx] -= d_m[d_idx]*ay + s_fz[s_idx] -= d_m[d_idx]*az class AkinciRigidFluidCoupling(Equation): diff --git a/pysph/sph/scheme.py b/pysph/sph/scheme.py index 0355d81d0..3e288d27c 100644 --- a/pysph/sph/scheme.py +++ b/pysph/sph/scheme.py @@ -355,7 +355,7 @@ def consume_user_options(self, options): self.configure(**data) def get_timestep(self, cfl=0.5): - return cfl*self.h0/self.c0 + return cfl * self.h0 / self.c0 def configure_solver(self, kernel=None, integrator_cls=None, extra_steppers=None, **kw): @@ -560,9 +560,9 @@ def consume_user_options(self, options): self.configure(**data) def get_timestep(self, cfl=0.25): - dt_cfl = cfl * self.h0/self.c0 + dt_cfl = cfl * self.h0 / self.c0 if self.nu > 1e-12: - dt_viscous = 0.125 * self.h0**2/self.nu + dt_viscous = 0.125 * self.h0**2 / self.nu else: dt_viscous = 1.0 dt_force = 1.0 @@ -702,6 +702,273 @@ def setup_properties(self, particles, clean=True): pa.set_output_arrays(output_props) +class BeadChainScheme(Scheme): + def __init__(self, fluids, solids, fibers, dim, k=0.0, lim=0.5, + tag=100, gx=0.0, gy=0.0, gz=0.0, alpha=0.0, tdamp=0.0, + direct=True): + self.fluids = fluids + self.solids = solids + self.fibers = fibers + self.dim = dim + self.solver = None + self.rho0 = None + self.c0 = None + self.pb = None + self.p0 = None + self.nu = None + self.h0 = None + self.dx = None + self.E = None + self.A = None + self.Ip = None + self.J = None + self.d = None + self.lim = lim + self.k = k + self.tag = tag + self.gx = gx + self.gy = gy + self.gz = gz + self.alpha = alpha + self.tdamp = tdamp + self.direct = direct + + def add_user_options(self, group): + group.add_argument( + "--tdamp", action="store", type=float, dest="tdamp", + default=self.tdamp, + help="Time for which the accelerations are damped." + ) + + def consume_user_options(self, options): + vars = ['tdamp'] + data = dict((var, getattr(options, var)) for var in vars) + self.configure(**data) + + def get_timestep(self): + from math import sqrt + dt_cfl = 0.4 * self.h0 / (1.1 * self.c0) + dt_viscous = 0.125 * self.h0**2 / self.nu + g = sqrt(self.gx**2 + self.gy**2 + self.gz**2) + dt_force = 0.25 * sqrt(self.h0 / (g + 0.001)) + dt_tension = 0.5 * self.dx * sqrt(self.rho0 / self.E) + dt_bending = 0.5 * self.dx**2 * sqrt(self.rho0 * self.A / + (self.E * 2 * self.Ip)) + print("dt_cfl: %g" % dt_cfl) + print("dt_viscous: %g" % dt_viscous) + print("dt_force: %g" % dt_force) + print("dt_tension: %g" % dt_tension) + print("dt_bending: %g" % dt_bending) + + fiber_dt = min(dt_tension, dt_bending) + dt = min(dt_cfl, dt_viscous, dt_force) + + if self.direct: + dt = min(dt_cfl, dt_viscous, dt_force, dt_tension, dt_bending) + if (dt / fiber_dt) > 10: + print("Time step ratio is %g." % (dt / fiber_dt)) + print("You might want to speed up the simulation using an") + print("inner loop. Set it up as") + print("def create_tools(self):") + print("return [FiberIntegrator(self.particles,") + print(" self.scheme,") + print(" self.domain,") + print(" D=)]") + else: + dt = min(dt_cfl, dt_viscous, dt_force) + print("Time step ratio is %g" % (dt / fiber_dt)) + return [dt, fiber_dt] + + def configure_solver(self, kernel=None, integrator_cls=None, + extra_steppers=None, **kw): + """Configure the solver to be generated. + + Parameters + ---------- + kernel : Kernel instance. + Kernel to use, if none is passed a default one is used. + integrator_cls : pysph.sph.integrator.Integrator + Integrator class to use, use sensible default if none is + passed. + extra_steppers : dict + Additional integration stepper instances as a dict. + **kw : extra arguments + Any additional keyword args are passed to the solver instance. + + """ + from pysph.base.kernels import CubicSpline + from pysph.sph.integrator_step import TransportVelocityStep + from pysph.sph.integrator import EPECIntegrator + if kernel is None: + kernel = CubicSpline(dim=self.dim) + steppers = {} + if extra_steppers is not None: + steppers.update(extra_steppers) + + step_cls = TransportVelocityStep + for f in self.fluids + self.fibers: + if f not in steppers: + steppers[f] = step_cls() + + cls = integrator_cls if integrator_cls is not None else EPECIntegrator + integrator = cls(**steppers) + + [self.dt, self.fiber_dt] = self.get_timestep() + + from pysph.solver.solver import Solver + self.solver = Solver(dim=self.dim, integrator=integrator, + kernel=kernel, dt=self.dt, **kw) + + def get_equations(self): + from pysph.sph.equation import Group + from pysph.sph.wc.transport_velocity import ( + SummationDensity, StateEquation, MomentumEquationPressureGradient, + MomentumEquationViscosity, MomentumEquationArtificialStress, + SolidWallPressureBC, SolidWallNoSlipBC, SetWallVelocity, + FiberViscousTraction) + from pysph.sph.fiber.utils import ( + HoldPoints, VelocityGradient, ComputeDistance, Contact) + from pysph.sph.fiber.beadchain import (EBGVelocityReset, Friction, + Tension, Bending) + + equations = [] + all = self.fluids + self.solids + self.fibers + + g1 = [] + for fluid in self.fluids: + g1.append(SummationDensity(dest=fluid, sources=all)) + + for fiber in self.fibers: + g1.append(EBGVelocityReset(dest=fiber, sources=None)) + equations.append(Group(equations=g1, real=False)) + + g2 = [] + for fluid in self.fluids: + g2.append(StateEquation( + dest=fluid, sources=None, p0=self.p0, rho0=self.rho0, b=1.0)) + + for solid in self.solids: + g2.append(SetWallVelocity(dest=solid, sources=self.fluids)) + + for fiber in self.fibers: + g2.append(SetWallVelocity(dest=fiber, sources=self.fluids)) + + equations.append(Group(equations=g2, real=False)) + + g3 = [] + for solid in self.solids: + g3.append(SolidWallPressureBC( + dest=solid, sources=self.fluids, b=1.0, + rho0=self.rho0, p0=self.p0, gx=self.gx, gy=self.gy, + gz=self.gz)) + for fiber in self.fibers: + g3.append(SolidWallPressureBC( + dest=fiber, sources=self.fluids, b=1.0, + rho0=self.rho0, p0=self.p0, gx=self.gx, gy=self.gy, + gz=self.gz)) + + equations.append(Group(equations=g3, real=False)) + + g4 = [] + for fiber in self.fibers: + g4.append(ComputeDistance(dest=fiber, sources=[fiber])) + g4.append(VelocityGradient(dest=fiber, sources=all)) + + equations.append(Group(equations=g4)) + + g5 = [] + for fluid in self.fluids: + g5.append(MomentumEquationPressureGradient( + dest=fluid, sources=all, pb=self.pb, gx=self.gx, gy=self.gy, + gz=self.gz, tdamp=self.tdamp)) + if self.nu > 0.0: + g5.append(MomentumEquationViscosity( + dest=fluid, sources=self.fluids, nu=self.nu)) + g5.append(SolidWallNoSlipBC( + dest=fluid, sources=self.solids + self.fibers, + nu=self.nu)) + g5.append(MomentumEquationArtificialStress( + dest=fluid, sources=self.fluids)) + + for fiber in self.fibers: + g5.append(MomentumEquationPressureGradient( + dest=fiber, sources=all, pb=0.0, gx=self.gx, gy=self.gy, + gz=self.gz, tdamp=self.tdamp)) + if self.nu > 0.0: + g5.append(FiberViscousTraction( + dest=fiber, sources=self.fluids, nu=self.nu)) + if self.direct: + g5.append(Tension(dest=fiber, + sources=None, + ea=self.E*self.A)) + g5.append(Bending(dest=fiber, + sources=None, + ei=self.E*self.Ip)) + g5.append(Contact(dest=fiber, + sources=self.fibers, + E=self.E, + d=self.d, + dim=self.dim, + k=self.k, + lim=self.lim, + eta0=self.rho0*self.nu, + dt=self.dt)) + + g5.append(Friction( + dest=fiber, sources=None, J=self.J, dx=self.dx, + mu=self.nu * self.rho0, d=self.d)) + + equations.append(Group(equations=g5)) + + g6 = [] + for fiber in self.fibers: + g6.append(HoldPoints(dest=fiber, sources=None, tag=100)) + + equations.append(Group(equations=g6)) + + return equations + + def setup_properties(self, particles, clean=True): + particle_arrays = dict([(p.name, p) for p in particles]) + + props = ('V', 'wf', 'uf', 'vf', 'wg', 'wij', 'vg', 'ug', + 'awhat', 'avhat', 'auhat', 'vhat', 'what', 'uhat', 'vmag2', + 'arho', 'rho0', 'holdtag') + fprops = ('lprev', 'lnext', + 'phi0', 'fidx', 'phifrac', 'phi0', 'fractag', 'holdtag', + 'rxnext', 'rynext', 'rznext', 'rnext', + 'rxprev', 'ryprev', 'rzprev', 'rprev', + 'Fx', 'Fy', 'Fz', + 'eu', 'ev', 'ew', + 'ex', 'ey', 'ez' + 'dudx', 'dudy', 'dudz', + 'dvdx', 'dvdy', 'dvdz', + 'dwdx', 'dwdy', 'dwdz') + output_props = ['x', 'y', 'z', 'u', 'v', 'w', 'rho', 'm', 'p', + 'pid', 'gid'] + + for fluid in self.fluids: + pa = particle_arrays[fluid] + for name in props: + fluid.add_property(name) + self._ensure_properties(pa, props, clean) + pa.set_output_arrays(output_props) + + for solid in self.solids: + pa = particle_arrays[solid] + for name in props: + solid.add_property(name) + self._ensure_properties(pa, props, clean) + pa.set_output_arrays(output_props) + + for fiber in self.fibers: + pa = particle_arrays[fiber] + for name in props + fprops: + fiber.add_property(name) + self._ensure_properties(pa, props, clean) + pa.set_output_arrays(output_props) + + class AdamiHuAdamsScheme(TVFScheme): """This is a scheme similiar to that in the paper: @@ -742,7 +1009,7 @@ def add_user_options(self, group): ) def attributes_changed(self): - self.B = self.c0*self.c0*self.rho0/self.gamma + self.B = self.c0 * self.c0 * self.rho0 / self.gamma def consume_user_options(self, options): vars = ['alpha', 'tdamp', 'gamma'] @@ -1494,7 +1761,7 @@ def get_equations(self): g4 = [] for fluid in self.fluids: - g4.append(SummationDensity(fluid, self.fluids+self.solids)) + g4.append(SummationDensity(fluid, self.fluids + self.solids)) equations.append(Group(g4)) g5 = [] @@ -1503,7 +1770,7 @@ def get_equations(self): equations.append(Group(equations=g5)) g6 = [] - for elem in self.fluids+self.solids: + for elem in self.fluids + self.solids: g6.append(IdealGasEOS(elem, sources=None, gamma=self.gamma)) equations.append(Group(equations=g6)) @@ -1574,10 +1841,10 @@ def setup_properties(self, particles, clean=True): import numpy particle_arrays = dict([(p.name, p) for p in particles]) required_props = [ - 'x', 'y', 'z', 'u', 'v', 'w', 'rho', 'h', 'm', 'cs', 'p', - 'e', 'au', 'av', 'aw', 'arho', 'ae', 'am', 'ah', 'x0', 'y0', - 'z0', 'u0', 'v0', 'w0', 'rho0', 'e0', 'h0', 'div', 'h0', - 'wij', 'htmp', 'logrho'] + 'x', 'y', 'z', 'u', 'v', 'w', 'rho', 'h', 'm', 'cs', 'p', + 'e', 'au', 'av', 'aw', 'arho', 'ae', 'am', 'ah', 'x0', 'y0', + 'z0', 'u0', 'v0', 'w0', 'rho0', 'e0', 'h0', 'div', 'h0', + 'wij', 'htmp', 'logrho'] dummy = get_particle_array(additional_props=required_props, name='junk') diff --git a/pysph/sph/wc/transport_velocity.py b/pysph/sph/wc/transport_velocity.py index 66db2d8b0..208c4a0d7 100644 --- a/pysph/sph/wc/transport_velocity.py +++ b/pysph/sph/wc/transport_velocity.py @@ -15,7 +15,7 @@ """ from pysph.sph.equation import Equation -from math import sin, pi +from math import pi, sin # constants M_PI = pi @@ -638,6 +638,54 @@ def loop(self, d_idx, s_idx, d_m, d_rho, s_rho, d_V, s_V, d_aw[d_idx] += tmp * (d_w[d_idx] - s_wg[s_idx]) +class FiberViscousTraction(Equation): + r"""**Fiber boundary condition** + """ + + def __init__(self, dest, sources, nu): + r""" + Parameters + ---------- + nu : float + kinematic viscosity + """ + + self.nu = nu + super(FiberViscousTraction, self).__init__(dest, sources) + + def initialize(self, d_idx, d_au, d_av, d_aw): + d_au[d_idx] = 0.0 + d_av[d_idx] = 0.0 + d_aw[d_idx] = 0.0 + + def loop(self, d_idx, s_idx, d_m, s_m, d_rho, s_rho, d_V, s_V, + s_u, s_v, s_w, + d_au, d_av, d_aw, + d_ug, d_vg, d_wg, + DWIJ, R2IJ, EPS, XIJ): + + # averaged shear viscosity Eq. (6). + etai = self.nu * d_rho[d_idx] + etaj = self.nu * s_rho[s_idx] + + etaij = 2 * (etai * etaj)/(etai + etaj) + + # particle volumes; d_V is inverse volume + Vi = 1./d_V[d_idx] + Vj = 1./s_V[s_idx] + Vi2 = Vi * Vi + Vj2 = Vj * Vj + + # scalar part of the kernel gradient + Fij = XIJ[0]*DWIJ[0] + XIJ[1]*DWIJ[1] + XIJ[2]*DWIJ[2] + + tmp = 1./d_m[d_idx] * (Vi2 + Vj2) * (etaij * Fij/(R2IJ + EPS)) + + d_au[d_idx] += tmp * (d_ug[d_idx]-s_u[s_idx]) + d_av[d_idx] += tmp * (d_vg[d_idx]-s_v[s_idx]) + d_aw[d_idx] += tmp * (d_wg[d_idx]-s_w[s_idx]) + + class SolidWallPressureBC(Equation): r"""**Solid wall pressure boundary condition** [Adami2012]_