tadd BC parameter to file/data structures, add utility functions - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
 (HTM) git clone git://src.adamsgaard.dk/sphere
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
 (DIR) commit b456fd10ab9a15540b915c36ee2d3ed08ca44faf
 (DIR) parent 6b42a6886653c7e1e620ca470ffa953820ca6de0
 (HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
       Date:   Thu, 11 Aug 2016 16:01:06 -0700
       
       add BC parameter to file/data structures, add utility functions
       
       Diffstat:
         M python/sphere.py                    |     131 ++++++++++++++++++++++++++++---
         M src/datatypes.h                     |       4 ++++
         M src/file_io.cpp                     |      10 ++++++++++
       
       3 files changed, 133 insertions(+), 12 deletions(-)
       ---
 (DIR) diff --git a/python/sphere.py b/python/sphere.py
       t@@ -329,11 +329,21 @@ class sim:
                    self.p_mod_f = numpy.zeros(1, dtype=numpy.float64)  # Frequency [Hz]
                    self.p_mod_phi = numpy.zeros(1, dtype=numpy.float64) # Shift [rad]
        
       +            if self.cfd_solver[0] == 1:  # Darcy solver
       +                # Boundary conditions at the sides of the fluid grid
       +                # 0: Dirichlet
       +                # 1: Neumann
       +                # 2: Periodic (default)
       +                self.bc_xn = numpy.ones(1, dtype=numpy.int32)*3  # Neg. x bc
       +                self.bc_xp = numpy.ones(1, dtype=numpy.int32)*3  # Pos. x bc
       +                self.bc_yn = numpy.ones(1, dtype=numpy.int32)*3  # Neg. y bc
       +                self.bc_yp = numpy.ones(1, dtype=numpy.int32)*3  # Pos. y bc
       +
                    # Boundary conditions at the top and bottom of the fluid grid
                    # 0: Dirichlet (default)
                    # 1: Neumann free slip
       -            # 2: Neumann no slip
       -            # 3: Periodic
       +            # 2: Neumann no slip (Navier Stokes), Periodic (Darcy)
       +            # 3: Periodic (Navier-Stokes solver only)
                    # 4: Constant flux (Darcy solver only)
                    self.bc_bot = numpy.zeros(1, dtype=numpy.int32)
                    self.bc_top = numpy.zeros(1, dtype=numpy.int32)
       t@@ -345,16 +355,6 @@ class sim:
                    self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
                    self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
        
       -            # Boundary conditions at the top and bottom of the fluid grid
       -            # 0: Dirichlet
       -            # 1: Neumann
       -            # 2: Periodic (default)
       -            # 3: Constant flux (Darcy solver only)
       -            self.bc_xn = numpy.zeros(1, dtype=numpy.int32)*2  # Neg. x boundary
       -            self.bc_xp = numpy.zeros(1, dtype=numpy.int32)*2  # Pos. x boundary
       -            self.bc_yn = numpy.zeros(1, dtype=numpy.int32)*2  # Neg. y boundary
       -            self.bc_yp = numpy.zeros(1, dtype=numpy.int32)*2  # Pos. y boundary
       -
                    ## Solver parameters
        
                    # Navier-Stokes
       t@@ -739,6 +739,18 @@ class sim:
                        elif self.k_c != other.k_c:
                            print(88)
                            return(88)
       +                elif (self.bc_xn != other.bc_xn):
       +                    print(92)
       +                    return 92
       +                elif (self.bc_xp != other.bc_xp):
       +                    print(93)
       +                    return 93
       +                elif (self.bc_yn != other.bc_yn):
       +                    print(94)
       +                    return 94
       +                elif (self.bc_yp != other.bc_yp):
       +                    print(95)
       +                    return 95
        
                if (self.color != other.color).any():
                    print(90)
       t@@ -1155,6 +1167,16 @@ class sim:
                            self.p_mod_phi =\
                                    numpy.fromfile(fh, dtype=numpy.float64, count=1)
        
       +                    if self.version >= 2.12 and self.cfd_solver == 1:
       +                        self.bc_xn =\
       +                            numpy.fromfile(fh, dtype=numpy.int32, count=1)
       +                        self.bc_xp =\
       +                            numpy.fromfile(fh, dtype=numpy.int32, count=1)
       +                        self.bc_yn =\
       +                            numpy.fromfile(fh, dtype=numpy.int32, count=1)
       +                        self.bc_yp =\
       +                            numpy.fromfile(fh, dtype=numpy.int32, count=1)
       +
                            self.bc_bot =\
                                    numpy.fromfile(fh, dtype=numpy.int32, count=1)
                            self.bc_top =\
       t@@ -1405,6 +1427,12 @@ class sim:
                        fh.write(self.p_mod_f.astype(numpy.float64))
                        fh.write(self.p_mod_phi.astype(numpy.float64))
        
       +                if self.cfd_solve[0] == 1:  # Sides only adjustable with Darcy
       +                    fh.write(self.bc_xn.astype(numpy.int32))
       +                    fh.write(self.bc_xp.astype(numpy.int32))
       +                    fh.write(self.bc_yn.astype(numpy.int32))
       +                    fh.write(self.bc_yp.astype(numpy.int32))
       +
                        fh.write(self.bc_bot.astype(numpy.int32))
                        fh.write(self.bc_top.astype(numpy.int32))
                        fh.write(self.free_slip_bot.astype(numpy.int32))
       t@@ -3358,6 +3386,11 @@ class sim:
                    self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
                    self.k_c = numpy.ones(1, dtype=numpy.float64)*4.6e-10
        
       +            self.bc_xn = numpy.ones(1, dtype=numpy.int32)*3
       +            self.bc_xp = numpy.ones(1, dtype=numpy.int32)*3
       +            self.bc_yn = numpy.ones(1, dtype=numpy.int32)*3
       +            self.bc_yp = numpy.ones(1, dtype=numpy.int32)*3
       +
                else:
                    raise Exception('Value of cfd_solver not understood (' + \
                            str(self.cfd_solver[0]) + ')')
       t@@ -3465,6 +3498,80 @@ class sim:
                self.bc_top[0] = 4
                self.bc_top_flux[0] = specific_flux
        
       +    def setFluidXFixedPressure(self):
       +        '''
       +        Set the X boundaries of the fluid domain to follow the fixed pressure
       +        value (Dirichlet) boundary condition.
       +
       +        This is not the default behavior for the boundary. See also
       +        :func:`setFluidXFixedPressure()`,
       +        :func:`setFluidXNoFlow()`, and
       +        :func:`setFluidXPeriodic()` (default)
       +        '''
       +        self.bc_xn[0] = 0
       +        self.bc_xp[0] = 0
       +
       +    def setFluidXNoFlow(self):
       +        '''
       +        Set the X boundaries of the fluid domain to follow the no-flow
       +        (Neumann) boundary condition.
       +
       +        This is not the default behavior for the boundary. See also
       +        :func:`setFluidXFixedPressure()`,
       +        :func:`setFluidXNoFlow()`, and
       +        :func:`setFluidXPeriodic()` (default)
       +        '''
       +        self.bc_xn[0] = 1
       +        self.bc_xp[0] = 1
       +
       +    def setFluidXPeriodic(self):
       +        '''
       +        Set the X boundaries of the fluid domain to follow the periodic
       +        (cyclic) boundary condition.
       +
       +        This is the default behavior for the boundary. See also
       +        :func:`setFluidXFixedPressure()` and
       +        :func:`setFluidXNoFlow()`
       +        '''
       +        self.bc_xn[0] = 2
       +        self.bc_xp[0] = 2
       +
       +    def setFluidYFixedPressure(self):
       +        '''
       +        Set the Y boundaries of the fluid domain to follow the fixed pressure
       +        value (Dirichlet) boundary condition.
       +
       +        This is not the default behavior for the boundary. See also
       +        :func:`setFluidYNoFlow()` and
       +        :func:`setFluidYPeriodic()` (default)
       +        '''
       +        self.bc_yn[0] = 0
       +        self.bc_yp[0] = 0
       +
       +    def setFluidYNoFlow(self):
       +        '''
       +        Set the Y boundaries of the fluid domain to follow the no-flow
       +        (Neumann) boundary condition.
       +
       +        This is not the default behavior for the boundary. See also
       +        :func:`setFluidYFixedPressure()` and
       +        :func:`setFluidYPeriodic()` (default)
       +        '''
       +        self.bc_yn[0] = 1
       +        self.bc_yp[0] = 1
       +
       +    def setFluidYPeriodic(self):
       +        '''
       +        Set the Y boundaries of the fluid domain to follow the periodic
       +        (cyclic) boundary condition.
       +
       +        This is the default behavior for the boundary. See also
       +        :func:`setFluidYFixedPressure()` and
       +        :func:`setFluidYNoFlow()`
       +        '''
       +        self.bc_yn[0] = 2
       +        self.bc_yp[0] = 2
       +
            def setPermeabilityPrefactor(self, k_c, verbose=True):
                '''
                Set the permeability prefactor from Goren et al 2011, eq. 24. The
 (DIR) diff --git a/src/datatypes.h b/src/datatypes.h
       t@@ -164,6 +164,10 @@ struct Darcy {
            Float   p_mod_A;        // Pressure modulation amplitude at top
            Float   p_mod_f;        // Pressure modulation frequency at top
            Float   p_mod_phi;      // Pressure modulation phase at top
       +    int     bc_xn;          // 0: Dirichlet, 1: Neumann, 3: Periodic
       +    int     bc_xp;          // 0: Dirichlet, 1: Neumann, 3: Periodic
       +    int     bc_yn;          // 0: Dirichlet, 1: Neumann, 3: Periodic
       +    int     bc_yp;          // 0: Dirichlet, 1: Neumann, 3: Periodic
            int     bc_bot;         // 0: Dirichlet, 1: Neumann
            int     bc_top;         // 0: Dirichlet, 1: Neumann
            int     free_slip_bot;  // 0: no, 1: yes
 (DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
       t@@ -294,6 +294,7 @@ void DEM::readbin(const char *target)
                    ifs.read(as_bytes(ns.p_mod_f), sizeof(Float));
                    ifs.read(as_bytes(ns.p_mod_phi), sizeof(Float));
        
       +            ifs.read(as_bytes(ns.bc_top), sizeof(int));
                    ifs.read(as_bytes(ns.bc_bot), sizeof(int));
                    ifs.read(as_bytes(ns.bc_top), sizeof(int));
                    ifs.read(as_bytes(ns.free_slip_bot), sizeof(int));
       t@@ -364,6 +365,10 @@ void DEM::readbin(const char *target)
                    ifs.read(as_bytes(darcy.p_mod_f), sizeof(Float));
                    ifs.read(as_bytes(darcy.p_mod_phi), sizeof(Float));
        
       +            ifs.read(as_bytes(darcy.bc_xn), sizeof(int));
       +            ifs.read(as_bytes(darcy.bc_xp), sizeof(int));
       +            ifs.read(as_bytes(darcy.bc_yn), sizeof(int));
       +            ifs.read(as_bytes(darcy.bc_yp), sizeof(int));
                    ifs.read(as_bytes(darcy.bc_bot), sizeof(int));
                    ifs.read(as_bytes(darcy.bc_top), sizeof(int));
                    ifs.read(as_bytes(darcy.free_slip_bot), sizeof(int));
       t@@ -659,6 +664,11 @@ void DEM::writebin(const char *target)
                        ofs.write(as_bytes(darcy.p_mod_f), sizeof(Float));
                        ofs.write(as_bytes(darcy.p_mod_phi), sizeof(Float));
        
       +                ofs.write(as_bytes(darcy.bc_xn), sizeof(int));
       +                ofs.write(as_bytes(darcy.bc_xp), sizeof(int));
       +                ofs.write(as_bytes(darcy.bc_yn), sizeof(int));
       +                ofs.write(as_bytes(darcy.bc_yp), sizeof(int));
       +                ofs.write(as_bytes(darcy.bc_top), sizeof(int));
                        ofs.write(as_bytes(darcy.bc_bot), sizeof(int));
                        ofs.write(as_bytes(darcy.bc_top), sizeof(int));
                        ofs.write(as_bytes(darcy.free_slip_bot), sizeof(int));