Merge pull request #48 from master-csmi/35-implement-the-sath-scheme

update

Library/SATh.py

@@ -24,6 +24,9 @@ def __init__(self, **kwargs):
         "Theta_max":1,
         "Scheme": "SATh", #"Theta"
         "Theta_computation":"Fixed-Point", #"Newton"
+        "Theta_choice_epsilon":1e-6,
+        "Method_convergence_epsilon":1e-6,
+        "Timer":False,
         "Path": "pictures",
         "Boundary":"constant"
         }
@@ -54,6 +57,9 @@ def __init__(self, **kwargs):
         else:
             print("Wrong Scheme type in the inputs")
 
+    def print_params(self):
+        return self.default_dict
+
     def sol_ex(self):
         return self.funcs.exact(self.mesh.nodes, self.params, self.tf)
 
@@ -94,6 +100,11 @@ def interpol_quadrature(self):  #use with tau=1/2 for standard linear interpolat
         #return ((self.env.mesh.dt/2)*self.u_down + (1 - self.env.mesh.dt/2)*self.u_up)/self.env.mesh.dt
         return self.u_down/2 + self.u_up/2
 
+    def timer(self, keyword=None):
+        if keyword == "set":
+            print(int(self.env.tf%self.env.mesh.dt) * "-")
+        else:
+            print("-",end="")
 
     def theta_choice(self, v, w, epsilon=1e-6):
         if np.abs(w) > epsilon:
@@ -135,26 +146,29 @@ def Newton(self, epsilon=1e-6):
         while np.linalg.norm(w - w_) > epsilon:
             for i in range(1,self.env.mesh.Nx +1):
                 w_[i] = w[i]
+
+                self.thetas[i] = self.theta_choice(v[i],w[i],epsilon=self.env.default_dict["Theta_choice_epsilon"])
+
                 #compute step (linear system resolution)
-                s = np.linalg.solve(self.compute_J(w,v,i), -self.F(w,v,i))
+                s = np.linalg.solve(self.compute_J(w,v,i, epsilon=self.env.default_dict["Theta_choice_epsilon"]), -self.F(w,v,i))
                 #iterate Newton problem
                 w[i] += s[0]
                 v[i] += s[1]
 
         return [w,v]
-    
+
 
     def fixed_point_iter(self, theta_left, u_down_left, u_down,
                                       w_left, lam, epsilon=1e-6):
         w = w_left #
         w_ = epsilon + w  #Just to initialize
-        theta = self.theta_st
+        theta = self.theta_st   #Try wit hthe previous time step theta
 
-        while np.abs(w_ - w) > epsilon:
+        while np.abs(w_ - w) >= epsilon:
             w_ = w
             w = (u_down - theta_left*w_left - (u_down_left + w_left)) / (1 + lam*theta)
             v = (-lam/2) * (u_down - theta_left**2 * w_left - (u_down_left + w_left) + theta**2*w)
-            theta = self.theta_choice(v,w)
+            theta = self.theta_choice(v,w, epsilon=self.env.default_dict["Theta_choice_epsilon"])
 
         return theta, w
 
@@ -168,13 +182,13 @@ def update_thetas(self):
 
             for i in range(1, self.env.mesh.Nx +1):
                 theta, w_left = self.fixed_point_iter(self.thetas[i-1], self.u_down[i-1],
-                                    self.u_down[i], w_left, self.lam)
+                                    self.u_down[i], w_left, self.lam, epsilon=self.env.default_dict["Method_convergence_epsilon"])
                 self.thetas[i] = theta
 
         elif self.theta_computation == "Newton":
-            [w,v] = self.Newton()
+            [w,v] = self.Newton(epsilon=self.env.default_dict["Method_convergence_epsilon"])
             for i in range(w.size):
-                self.thetas[i] = self.theta_choice(v[i],w[i])
+                self.thetas[i] = self.theta_choice(v[i],w[i], epsilon=self.env.default_dict["Theta_choice_epsilon"])
 
             self.u_til = v
 
@@ -186,22 +200,30 @@ def SATh_Scheme(self):
         self.u_up = np.empty_like(self.u_down)
         self.u_up[0] = self.u_down[0]
 
+        if self.env.default_dict["Timer"]==True:
+            self.timer("set")
+
         while (t<self.env.tf):
 
             if t!=0:
                 self.update_thetas()
+            #print(min(self.thetas), max(self.thetas))
 
             coefs = self.env.mesh.dt*(np.eye(self.env.mesh.Nx+1)-np.diag(self.thetas))
             A = self.env.mats.Iter_Mat(self.env.mesh, self.thetas, self.env.alpha, adaptive=True)
+            #A = sparse.linalg.LinearOperator((self.env.mesh.Nx+1,self.env.mesh.Nx+1), matvec=self.env.funcs.Iter_Func)
             b = self.u_down - self.env.alpha * coefs@(self.env.mats.Dx @ self.u_down)
             b[0] = self.u_down[0]
 
             self.u_up, _ = gmres(A, b)
+            #self.u_up = A.matvec(b)
             self.u_down = self.u_up.copy()
             self.u_til += self.u_up    #Only needed in the Newton method to update v (-> see line 176)
 
             t += self.env.mesh.dt
-            #print(t)
+            if self.env.default_dict["Timer"]==True:
+                self.timer()
+                #print(t)
 
         return self.u_up
 

Library/utilities.py

@@ -73,7 +73,14 @@ def Iter_Mat(self, mesh, theta, alpha, adaptive):
             A = Id + ((Dx * theta[:,np.newaxis]) * mesh.dt * alpha)
             A[0,0] = 1
             return sparse.csr_matrix(A)
+
+    def Iter_Func(self, mesh, theta, alpha, v):   #Function version of Iter_Mat in order to build a linear operator.
+                                               #For SATh, as in the case of the simple theta scheme the matrix is only needed to be built once for all.
+        for i in range(1,v.size):
+            v[i] = v[i]*(1+ alpha*theta[i]*mesh.dt/mesh.dx) - v[i-1] * alpha*theta[i]*mesh.dt/mesh.dx
 
+        return v
+
 
 class Functions():
     def __init__(self, mesh, params, tf, init_type):

args.json

@@ -12,6 +12,9 @@
     "Theta_max":1,
     "Scheme": "SATh",
     "Theta_computation": "Fixed-Point",
+    "Theta_choice_epsilon":1e-6,
+    "Method_convergence_epsilon":1e-6,
+    "Timer":0,
     "Path": "pictures",
     "Boundary":"constant"
 }

texdir/report.tex

@@ -233,12 +233,14 @@ \subsection{Self-Adaptive Theta Scheme for the advection equation}
 
 \begin{align*}
     \theta_i^{n+1} = \begin{cases}
-        max(\theta_{min}, \frac{\tilde{u}_i^{n+1} - u_i^n}{u_i^{n+1} - u_i^n} ) & \text{if } |u_i^{n+1} - u_i^n| > \epsilon \tag{4} \\
+        min(max(\theta_{min}, \frac{\tilde{u}_i^{n+1} - u_i^n}{u_i^{n+1} - u_i^n} ), 1) & \text{if } |u_i^{n+1} - u_i^n| > \epsilon \tag{4} \\
         \theta^* & \text{else} 
     \end{cases}
 \end{align*}
+
 We introduce $\theta_{min}$ and $\theta^*$ as new parameters of the numerical method.\\
 The value $\tilde{u}_i^{n+1}$ is an interpolation of the numerical solution between two time steps.\\
+
 \begin{figure}[H]
     \centering
     \includegraphics[width=0.7\textwidth]{indexmapu.png}
@@ -453,4 +455,11 @@ \section{Bibliography}
     \item Berzins and Furzeland, \textit{An adaptive theta method for the solution of stiff and nonstiff differential equations}, 1992
 \end{itemize}
 
+\begin{align*}
+    \theta_i^{n+1} = \begin{cases}
+        min(max(\theta_{min}, \frac{\tilde{u}_i^{n+1} - u_i^n}{u_i^{n+1} - u_i^n} ), 1) & \text{if } |u_i^{n+1} - u_i^n| > \epsilon \tag{4} \\
+        \theta^* & \text{else} 
+    \end{cases}
+\end{align*}
+
 \end{document}