tadded IO - ns2dfd - 2D finite difference Navier Stokes solver for fluid dynamics
 (HTM) git clone git://src.adamsgaard.dk/ns2dfd
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
 (DIR) commit 0a69461f671d81906ae29bc2bbed530b9aee906a
 (DIR) parent 02907fc07a66c16364379dbf5d92550083ce0730
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Sat,  1 Mar 2014 17:53:55 +0100
       
       added IO
       
       Diffstat:
         M ns2dfd.py                           |     155 ++++++++++++++++++++++++++-----
       
       1 file changed, 132 insertions(+), 23 deletions(-)
       ---
 (DIR) diff --git a/ns2dfd.py b/ns2dfd.py
       t@@ -4,17 +4,18 @@ import numpy
        
        class fluid:
        
       -    def __init__(sim = 'unnamed'):
       +    def __init__(name = 'unnamed'):
                '''
                A Navier-Stokes two-dimensional fluid flow simulation object. Most
                simulation values are assigned default values upon initialization.
       -        :param sim: Simulation identifier
       -        :type sim: str
       +        :param name: Simulation identifier
       +        :type name: str
                '''
       -        self.sim = sim
       +        self.sim_id = name
        
                init_grid()
                current_time()
       +        end_time()
                safety_factor()
                max_iterations()
                tolerance_criteria()
       t@@ -36,11 +37,13 @@ class fluid:
                :param dy: Grid cell height (meters)
                :type dy: float
                '''
       -        self.nx, self.ny = nx, ny
       -        self.dx, self.dy = dx, dy
       +        self.nx = numpy.asarray(nx)
       +        self.ny = numpy.asarray(ny)
       +        self.dx = numpy.asarray(dx)
       +        self.dy = numpy.asarray(dy)
       +        self.p = numpy.zeros((nx+1, ny+1))
                self.u = numpy.zeros((nx+2, ny+2))
                self.v = numpy.zeros((nx+2, ny+2))
       -        self.p = numpy.zeros((nx+1, ny+1))
        
            def current_time(t = 0.0):
                '''
       t@@ -48,16 +51,15 @@ class fluid:
                :param t: The current time value.
                :type t: float
                '''
       -        self.t_end = t_end
       +        self.t_end = numpy.asarray(t)
        
       -
       -    def end_time(t_end):
       +    def end_time(t_end = 1.0):
                '''
                Set the simulation end time.
                :param t_end: The time when to stop the simulation.
                :type t_end: float
                '''
       -        self.t_end = t_end
       +        self.t_end = numpy.asarray(t_end)
        
            def safety_factor(tau = 0.5):
                '''
       t@@ -66,7 +68,7 @@ class fluid:
                :param tau: Safety factor in ]0;1]
                :type tau: float
                '''
       -        self.tau = tau
       +        self.tau = numpy.asarray(tau)
                
            def max_iterations(itermax = 5000):
                '''
       t@@ -74,7 +76,7 @@ class fluid:
                :param itermax: Max. solution iterations in [1;inf[
                :type itermax: int
                '''
       -        self.itermax = itermax
       +        self.itermax = numpy.asarray(itermax)
        
            def tolerance_criteria(epsilon = 1.0e-4):
                '''
       t@@ -82,7 +84,7 @@ class fluid:
                :param epsilon: Criteria value
                :type epsilon: float
                '''
       -        self.epsilon = epsilon
       +        self.epsilon = numpy.asarray(epsilon)
        
            def relaxation_parameter(omega = 1.7):
                '''
       t@@ -92,7 +94,7 @@ class fluid:
                :param omega: Relaxation parameter value, in ]0;2[
                :type omega: float
                '''
       -        self.omega = omega
       +        self.omega = numpy.asarray(omega)
        
            def upwind_differencing_factor(gamma = 0.9):
                '''
       t@@ -101,7 +103,7 @@ class fluid:
                :param gamma: Upward differencing factor value, in ]0;1[
                :type gamma: float
                '''
       -        self.gamma = gamma
       +        self.gamma = numpy.asarray(gamma)
        
            def boundary_conditions(left = 1, right = 1, top = 1, bottom = 1):
                '''
       t@@ -110,10 +112,10 @@ class fluid:
                :param left, right, top, bottom: The wall to specify the BC for
                :type left, right, top, bottom: int
                '''
       -        self.w_left = left
       -        self.w_right = right
       -        self.w_top = top
       -        self.w_bottom = bottom
       +        self.w_left = numpy.asarray(left)
       +        self.w_right = numpy.asarray(right)
       +        self.w_top = numpy.asarray(top)
       +        self.w_bottom = numpy.asarray(bottom)
        
            def reynolds_number(re = 100):
                '''
       t@@ -121,7 +123,7 @@ class fluid:
                :param re: Reynolds number in ]0;infty[
                :type re: float
                '''
       -        self.re = re
       +        self.re = numpy.asarray(re)
        
            def gravity(gx = 0.0, gy = 0.0):
                '''
       t@@ -132,6 +134,113 @@ class fluid:
                    downward.
                :type gy: float
                '''
       -        self.gx, self.gy = gx, gy
       -
       +        self.gx = numpy.asarray(gx)
       +        self.gy = numpy.asarray(gy)
       +
       +    def read(path, verbose = True):
       +        '''
       +        Read data file from disk.
       +        :param path: Path to data file
       +        :type path: str
       +        '''
       +        fh = None
       +        try:
       +            targetfile = folder + '/' + self.sim_id
       +            if verbose == True:
       +                print('Input file: ' + targetfile)
       +            fh = open(targetfile, 'rb')
       +
       +            self.t     = numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +            self.t_end = numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +            self.tau   = numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +            
       +            self.itermax = numpy.fromfile(fh, dtype=numpy.int32, count=1)
       +            self.epsilon = numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +            self.omega   = numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +            self.gamma   = numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +
       +            self.gx = numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +            self.gy = numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +            self.re = numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +
       +            self.w_left   = numpy.fromfile(fh, dtype=numpy.int32, count=1)
       +            self.w_right  = numpy.fromfile(fh, dtype=numpy.int32, count=1)
       +            self.w_top    = numpy.fromfile(fh, dtype=numpy.int32, count=1)
       +            self.w_bottom = numpy.fromfile(fh, dtype=numpy.int32, count=1)
       +
       +            self.dx = numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +            self.dy = numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +            self.nx = numpy.fromfile(fh, dtype=numpy.int32, count=1)
       +            self.ny = numpy.fromfile(fh, dtype=numpy.int32, count=1)
       +
       +            init_grid(dx = self.dx, dy = self.dy, nx = self.nx, ny = self.ny)
       +            for j in range(self.ny+1):
       +                for i in range(self.nx+1):
       +                    self.p[i,j] = \
       +                            numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +
       +            for j in range(self.ny+2):
       +                for i in range(self.nx+2):
       +                    self.u[i,j] = \
       +                            numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +
       +            for j in range(self.ny+2):
       +                for i in range(self.nx+2):
       +                    self.v[i,j] = \
       +                            numpy.fromfile(fh, dtype=numpy.float64, count=1)
       +
       +        finally:
       +            if fh is not None:
       +                fh.close()
       +
       +    def write(verbose = True, folder = './'):
       +        '''
       +        Write the simulation parameters to disk so that the fluid flow solver
       +        can read it.
       +        '''
       +        fh = None
       +        try:
       +            targetfile = folder + '/' + self.sim_id
       +            if verbose == True:
       +                print('Output file: ' + targetfile)
       +            fh = open(targetfile, 'wb')
       +
       +            fh.write(self.t.astype(numpy.float64))
       +            fh.write(self.t_end.astype(numpy.float64))
       +            fh.write(self.tau.astype(numpy.float64))
       +            
       +            fh.write(self.itermax.astype(numpy.int32))
       +            fh.write(self.epsilon.astype(numpy.float64))
       +            fh.write(self.omega.astype(numpy.float64))
       +            fh.write(self.gamma.astype(numpy.float64))
       +
       +            fh.write(self.gx.astype(numpy.float64))
       +            fh.write(self.gy.astype(numpy.float64))
       +            fh.write(self.re.astype(numpy.float64))
       +
       +            fh.write(self.w_left.astype(numpy.int32))
       +            fh.write(self.w_right.astype(numpy.int32))
       +            fh.write(self.w_top.astype(numpy.int32))
       +            fh.write(self.w_bottom.astype(numpy.int32))
       +
       +            fh.write(self.dx.astype(numpy.float64))
       +            fh.write(self.dy.astype(numpy.float64))
       +            fh.write(self.nx.astype(numpy.int32))
       +            fh.write(self.ny.astype(numpy.int32))
       +
       +            for j in range(self.ny+1):
       +                for i in range(self.nx+1):
       +                    fh.write(self.p[i,j].astype(numpy.float64))
       +
       +            for j in range(self.ny+2):
       +                for i in range(self.nx+2):
       +                    fh.write(self.u[i,j].astype(numpy.float64))
       +
       +            for j in range(self.ny+2):
       +                for i in range(self.nx+2):
       +                    fh.write(self.v[i,j].astype(numpy.float64))
       +
       +        finally:
       +            if fh is not None:
       +                fh.close()