assignment11/assignment11.py at master · PGE323M/assignment11 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
import numpy as np
import yaml
import scipy.sparse
import scipy.sparse.linalg
import matplotlib.pyplot as plt

class OneDimReservoir():

    def __init__(self, inputs):
        '''
            Class for solving one-dimensional reservoir problems with
            finite differences.
        '''

        #stores input dictionary as class attribute, either read from a yaml file
        #or directly from a Python dictonary
        if isinstance(inputs, str):
            with open(inputs) as f:
                self.inputs = yaml.safe_load(f)
        else:
            self.inputs = inputs

        #computes delta_x
        self.Nx = self.inputs['numerical']['number of grids']['x']
        self.N = self.Nx
        if self.Nx < 2:
            raise ValueError('number of grids must be at least two')
        self.delta_x = self.inputs['reservoir']['length'] / float(self.Nx)

        #gets delta_t from inputs
        self.delta_t = self.inputs['numerical']['time step']

        #applies the initial reservoir pressues to self.p
        self.apply_initial_conditions()

        #calls fill matrix method (must be completely implemented to work)
        self.fill_matrices()

        #create an empty list for storing data if plots are requested
        if 'plots' in self.inputs:
            self.p_plot = []

        return

    def compute_transmissibility(self):
        '''
            Computes the transmissibility.
        '''

        raise NotImplementedError('Implement compute_transmissibility')


    def compute_accumulation(self):
        '''
            Computes the accumulation.
        '''

        raise NotImplementedError('Implement compute_accumulation')

    def fill_matrices(self):
        '''
            Fills the matrices T, B, and \vec{Q} and applies boundary
            conditions.
        '''

        raise NotImplementedError('Implement fill_matrices')


    def apply_boundary_conditions(self, T, Q):
        '''
            Applies boundary conditions (called by fill_matrices).
        '''

        N = self.Nx

        bc = self.inputs['boundary conditions']

        factor = self.inputs['conversion factor']
        for boundary in bc.values():
            if boundary['type'] not in ('prescribed pressure', 'prescribed flux'):
                raise ValueError('unsupported boundary type')
            if boundary['type'] == 'prescribed flux' and boundary['value'] != 0:
                raise ValueError('only zero prescribed flux is supported')

        if bc['left']['type'] == 'prescribed pressure':

            T[0, 0] = 3
            Q[0] = 2 * bc['left']['value']  * self.compute_transmissibility() * factor

        elif bc['left']['type'] == 'prescribed flux':

            T[0, 0] = 1
            Q[0] = 0

        if bc['right']['type'] == 'prescribed pressure':

            T[N-1, N-1] = 3
            Q[N-1] = 2 * bc['right']['value'] * self.compute_transmissibility() * factor

        elif bc['right']['type'] == 'prescribed flux':

            T[N-1, N-1] = 1
            Q[N-1] = 0

        return


    def apply_initial_conditions(self):
        '''
            Applies initial pressures to self.p
        '''

        N = self.N

        self.p = np.ones(N) * self.inputs['initial conditions']['pressure']

        return


    def solve_one_step(self):
        '''
            Solve one time step using either the implicit or explicit method
        '''

        raise NotImplementedError('Implement solve_one_step')


    def solve(self):
        '''
            Solves until "number of time steps"
        '''

        for i in range(self.inputs['numerical']['number of time steps']):
            self.solve_one_step()

            if 'plots' in self.inputs and i % self.inputs['plots']['frequency'] == 0:
                self.p_plot.append(self.get_solution().copy())

        return

    def plot(self):
        '''
           Crude plotting function.  Plots pressure as a function of grid block #
        '''

        if hasattr(self, 'p_plot'):
            for i in range(len(self.p_plot)):
                plt.plot(self.p_plot[i])

        return

    def get_solution(self):
        '''
            Returns solution vector
        '''
        return self.p.copy()