Implicit solver extracted to class

diff --git a/const.py b/const.py
index 20def04..e85d993 100644 (file)
--- a/const.py
+++ b/const.py
@@ -18,8 +18,8 @@ def get_characteristic(x):
         l = s.split('-')
         return int('-' + l[1])
 
-def timestep(step):
-    val = step / (C * sqrt(2))
+def timestep(val):
+    # val = step / (C * sqrt(2))
     c = get_characteristic(val)
     return power(10, double(c))
     

diff --git a/grid.py b/grid.py
index bd32b37..f63c0a4 100644 (file)
--- a/grid.py
+++ b/grid.py
@@ -10,11 +10,11 @@ class Grid(object):
     pass
 
 
-class Solver(object):
+class ExplicitSolver(object):
     '''Represents equation solver with time stepper    
     '''
 
-    def __init__(self, I, J, time, ht=None, ly=2e-5, lz=2e-5, lambda_=2e-6, C=3e8, scheme='explicit', solver='blitz', debug=False):
+    def __init__(self, I, J, time, ht=None, ly=2e-5, lz=2e-5, lambda_=2e-6, C=3e8, solver='blitz', debug=False):
 
         if solver == 'python':
             self._step = self._python_step
@@ -25,12 +25,7 @@ class Solver(object):
         else:
             raise ValueError('incorrect solver: %s' % solver)
 
-        if scheme == 'explicit':
-            self.solve = self._explicit_solve
-        elif scheme == 'implicit':
-            self.solve = self._implicit_solve
-        else:
-            raise ValueError('incorrect scheme type: %s' % scheme)
+        self.solve = self._explicit_solve
         
         self.C = C
         self.time = time
@@ -38,11 +33,7 @@ class Solver(object):
         self.hy = np.double(ly / I)
         self.hz = np.double(lz / J)
         if ht is None:
-            if scheme == 'explicit':
-                self.ht = timestep(min(self.hy, self.hz))
-            else:
-                # ???
-                self.ht = timestep(min(self.hy, self.hz))
+            self.ht = timestep(min(self.hy, self.hz) / (C * np.sqrt(2)))
         else:
             self.ht = ht
         self.K = int(self.time / self.ht)
@@ -122,6 +113,50 @@ class Solver(object):
             self._solved = True
         return self.u[2]
 
+
+class ImplicitSolver(object):
+    '''Represents equation solver with implicit difference scheme
+    '''
+    
+    def __init__(self, I, J, time, ht=None, ly=2e-5, lz=2e-5, lambda_=2e-6, C=3e8, debug=False):
+
+        self.solve = self._implicit_solve
+        
+        self.C = C
+        self.time = time
+        self.I, self.J, self.ly, self.lz, self.lambda_ = I, J, ly, lz, lambda_
+        self.hy = np.double(ly / I)
+        self.hz = np.double(lz / J)
+        if ht is None:
+            self.ht = timestep(min(self.hy, self.hz) / (C * np.sqrt(2)))
+        else:
+            self.ht = ht
+        self.K = int(self.time / self.ht)
+        self.u = np.zeros((3, I + 1, J + 1), 'd')
+
+        self._init_bc(0, 0)
+        self._init_bc(1, 1)
+        self._init_bc(2, 2)
+        self._solved = False
+        
+        if debug:
+            self._debug_print()
+
+    def _debug_print(self):
+        print "grid: %d x %d" % (self.I, self.J)
+        print "time:ht:K = %g:%g:%d" % (self.time, self.ht, self.K)
+        print
+        print "ly:hy:I = %g:%g:%d" % (self.ly, self.hy, self.I)
+        print "lz:hz:J = %g:%g:%d" % (self.lz, self.hz, self.J)
+        print
+        print "u size (in Kbytes):", ((self.u.size * self.u.itemsize) / 1024.)
+        print "extrapolated size (in Mbytes):", (float(self.I * self.J * self.K * self.u.itemsize) / (1024**2))
+        print
+
+    def _init_bc(self, k, time):
+        for i in xrange(self.u.shape[1]):
+            self.u[k, i, 0] = np.sin((np.pi * i) / self.I) * np.sin((2 * np.pi * self.C * self.ht * time) / self.lambda_)
+
     def _make_system(self):
         '''Формирует матрицу системы для прогонки (ее диагонали).
         '''
@@ -154,6 +189,5 @@ class Solver(object):
                     self.u[2, i, 1:-1] = sweep(a, b, c, f)
                 self.u[0] = self.u[1].copy()
                 self.u[1] = self.u[2].copy()
-                # self._init_bc(2, k + 1)
             self._solver = True
         return self.u[2]

diff --git a/main.py b/main.py
index 475701f..edb8a03 100755 (executable)
--- a/main.py
+++ b/main.py
@@ -64,14 +64,14 @@ def main():
     solveri = Solver(ny, nz, time, ht, scheme='implicit', debug=True)
     ui = solveri.solve()
 
-    solvere = Solver(ny, nz, time, ht, scheme='explicit')
-    ue = solvere.solve()
+    #solvere = Solver(ny, nz, time, ht, scheme='explicit')
+    #ue = solvere.solve()
     
     d = np.arange(0, solveri.J + 1)
     z = [solveri.hz * j for j in d]
     Ea = [Exm(const.LY/2, z_, time, 1000) for z_ in z]
     Ei = ui[ui.shape[0] / 2, :]
-    Ee = ue[ue.shape[0] / 2, :]
+    #Ee = ue[ue.shape[0] / 2, :]
 
 #     err = np.float64('0')
 #     for i in xrange(0, len(Ea)):
@@ -81,8 +81,8 @@ def main():
 
     a_line, = plt.plot(z, Ea, lw=1, color="red")
     i_line, = plt.plot(z, Ei, lw=1, color="green")
-    e_line, = plt.plot(z, Ee, lw=1, color="blue")
-    plt.legend((a_line, i_line, e_line), (u"Analitycal", u"Implicit scheme", u"Explicit scheme"))
+    #e_line, = plt.plot(z, Ee, lw=1, color="blue")
+    plt.legend((a_line, i_line), (u"Analitycal", u"Implicit scheme"))
     plt.grid(True)
     plt.xlim((min(z), max(z)))
     plt.show()