Skip to content

Commit

Permalink
more cleaning / advancement
Browse files Browse the repository at this point in the history
  • Loading branch information
@ARegardin
ARegardin committed Jun 18, 2024
1 parent 74cce09 commit faf8db3
Show file tree
Hide file tree
Showing 2 changed files with 133 additions and 113 deletions.
243 changes: 131 additions & 112 deletions Library/SATh.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,11 +9,12 @@
class Problem():
def __init__(self, **kwargs):

self.default_dict = {
self.params_dict = {
"a": 0,
"b": 1,
"Nx": 500,
"cfl": 5,
"f": "linear_advection",
"alpha": 1,
"mesh_type": "offset_constantDx",
"params": 0.2,
Expand All @@ -22,37 +23,39 @@ def __init__(self, **kwargs):
"Theta_st":.5,
"Theta_min":0.2,
"Theta_max":1,
"Scheme": "SATh", #"Theta"
"Scheme": "SATh", #/"Theta"
"Flux": "UP", #/"LF"
"Theta_computation":"Newton", #"Fixed-Point" has been abandoned
"Theta_choice_epsilon":1e-6,
"Theta_choice_method":"MinMax", #"Smoothing"
"Theta_choice_method":"MinMax", #/"Smoothing"
"Method_convergence_epsilon":1e-6,
"Timer":False,
"print_Newton_iter":False,
"Path": "pictures",
"Boundary":"constant"
"Boundary":"constant",
"tanh_factor":100,
}
#replacing by input values in the dict:
for key in kwargs.keys():
self.default_dict[key] = kwargs[key]
self.params_dict[key] = kwargs[key]

self.params = self.default_dict["params"]
self.tf = self.default_dict["tf"]
self.alpha = self.default_dict["alpha"]
self.params = self.params_dict["params"]
self.tf = self.params_dict["tf"]
self.alpha = self.params_dict["alpha"]

self.mesh = Mesh(self.default_dict["a"], self.default_dict["b"],
self.default_dict["Nx"], self.default_dict["cfl"],
self.default_dict["Boundary"], self.default_dict["mesh_type"])
self.mats = Matrices(self.mesh, self.default_dict["Boundary"], self.alpha)
self.funcs = Functions(self.mesh, self.params, self.tf, self.default_dict["init_func"])
self.mesh = Mesh(self.params_dict["a"], self.params_dict["b"],
self.params_dict["Nx"], self.params_dict["cfl"],
self.params_dict["Boundary"], self.params_dict["mesh_type"])
self.mats = Matrices(self.mesh, self.params_dict["Boundary"], self.alpha)
self.funcs = Functions(self.mesh, self.params, self.tf, self.params_dict["init_func"])

if self.default_dict["Scheme"] == "Theta": #Simple/Standard Theta Scheme
self.theta = self.default_dict["Theta_st"]
if self.params_dict["Scheme"] == "Theta": #Simple/Standard Theta Scheme
self.theta = self.params_dict["Theta_st"]
scheme = Theta_Scheme(self)
self.sol_num = scheme.solver()
self.sol_exc = self.sol_ex()

elif self.default_dict["Scheme"] == "SATh": #Self-Adaptive Theta Scheme
elif self.params_dict["Scheme"] == "SATh": #Self-Adaptive Theta Scheme
self.solver = SATh_Solver(self)
self.sol_num = self.solver.SATh_Scheme()
self.sol_exc = self.sol_ex()
Expand All @@ -61,156 +64,172 @@ def __init__(self, **kwargs):
print("Wrong Scheme type in the inputs")

def print_params(self):
return self.default_dict
return self.params_dict

def sol_ex(self):
return self.funcs.exact(self.mesh.nodes, self.params, self.tf)

def f(self, u):
if self.params_dict["f"] == "linear_advection":
return self.alpha * u



class SATh_Solver: #Works only for nonperiodical boundary conditions!
def __init__(self, env):
self.env = env
self.theta_st = env.default_dict["Theta_st"]
self.theta_min = env.default_dict["Theta_min"]
self.theta_st = env.params_dict["Theta_st"]
self.theta_min = env.params_dict["Theta_min"]
self.thetas = np.ones(env.mesh.Nx+1) * self.theta_st

self.u_down = self.env.funcs.init_sol.copy()
self.u_til = np.empty_like(self.u_down)
self.u_up = np.empty_like(self.u_down)
self.left_bound_val = self.env.funcs.init_sol[0]
self.u_up, self.u_til = np.empty_like(self.u_down), np.empty_like(self.u_down)
self.left_bound_val = self.env.funcs.init_sol[0] #!Change to self.bound_vals = [left,right]

self.theta_computation = self.env.default_dict["Theta_computation"]
self.w, self.v = np.empty_like(self.u_down), np.empty_like(self.u_down)

self.lam = env.mesh.dt/env.mesh.dx #CFL, here with the case of the flux function f(u) = au
#where a=1
#self.cyclicpermut_mat = sparse.csr_matrix(np.roll(np.eye(env.mesh.Nx), shift=-1, axis=0))
#self.d = lambda u: self.env.funcs.arrange_u(u,-1,bc=self.default_dict["Boundary"])

def timer(self, keyword=None):
if keyword == "set":
print(int(self.env.tf%self.env.mesh.dt) * "-")
self.theta_computation = self.env.params_dict["Theta_computation"]

if env.params_dict["Flux"]=="UP":
self.alpha = self.env.alpha
self.lam = env.mesh.dt/env.mesh.dx #CFL, here with the case of the flux function f(u) = au
elif env.params_dict["Flux"]=="LF":
self.alpha = self.alpha_LF()
self.lam = .5 * env.mesh.dt/env.mesh.dx
else: print("Unknown Flux type")


self.tanh_factor = self.env.params_dict["tanh_factor"]

def timer(self, action=None):
if action == "set":
print(int(self.env.tf//self.env.mesh.dt) * "_")
else:
print("-",end="")

def sech(self, x):
return 1/np.cosh(x)

def theta_choice(self, v, w, i, epsilon=1e-6):
if self.env.default_dict["Theta_choice_method"] == "MinMax":
if np.abs(w) > epsilon:
#print(np.abs(v/w > 1))
#if np.abs(v/w) > self.theta_min:
# print("theta=v/w",self.env.default_dict["Theta_computation"])
# print(max(self.theta_min, np.abs(v/w) ))
#return min(max(self.theta_min, np.abs(v/w) ), 1) #
self.thetas[i] = min(max(self.theta_min, np.abs(v/w) ), 1)
def alpha_LF(self):
pass

def theta_choice(self, i, epsilon=1e-6):
if self.env.params_dict["Theta_choice_method"] == "MinMax":
if np.abs(self.w[i]) > epsilon:
self.thetas[i] = min(max(self.theta_min, np.abs(self.v[i]/self.w[i]) ), 1)
else:
self.thetas[i] = self.theta_st

elif self.env.default_dict["Theta_choice_method"] == "Smoothing":
self.thetas[i] = .75 + .25 * np.tanh(100 * (-.5 + v/w))
elif self.env.params_dict["Theta_choice_method"] == "Smoothing":
if self.w[i]==0:
self.thetas[i] = 1
else:
self.thetas[i] = .75 + .25 * np.tanh(self.tanh_factor * (-.5 + self.v[i]/self.w[i]))

else :
print("Wrong Theta choice method type")


def F(self, w, v, i):
return np.array([
w[i] + self.lam * self.thetas[i] * w[i] + self.lam * (self.u_up[i] - self.thetas[i-1]*w[i-1] - self.u_up[i-1]),
v[i] + self.lam * .5 * self.thetas[i]**2 * w[i] + self.lam * .5 * (self.u_up[i] - self.thetas[i-1]**2 * w[i-1] - self.u_up[i-1])
])

def compute_J(self, w, v, i, epsilon=1e-6): #Compute the Jacobian
if self.env.default_dict["Theta_choice_method"] == "MinMax":
if np.abs(w[i])>epsilon :
if v[i]/w[i] > self.theta_min:
J = np.array([[1, self.lam], [-(v[i]**2)*self.lam / (2*(w[i]**2)), 1 + self.lam/w[i]]])
else :
J = np.array([[1+self.lam*self.theta_min, 0], [self.lam * (self.theta_min**2) / 2, 1]])
else :
J = np.array([[1+self.lam*self.theta_st, 0], [self.lam * (self.theta_st**2) / 2, 1]])

if self.env.default_dict["Theta_choice_method"] == "Smoothing":
J = np.array([[1 - 25 * v[i] * self.lam * np.sech(100*(-.5+ v[i]/w[i]))**2 / w[i] + self.lam * (3/4 + (1/4)*np.tanh(100*(-.5 + v[i]/w[i]))) ,
25 * self.lam * np.sech(100*(-.5 + v[i]/w[i]))**2],
[-(25 * v[i] * self.lam * np.sech(100*(-.5 + v[i]/w[i]))**2 * (3/4 + (1/4)*np.tanh(100*(-.5 + v[i]/w[i])))) / w + .5 * self.lam * (3/4 + (1/4)*np.tanh(100 * (-.5 + v[i]/w[i])))**2 ,
1 + 25 * self.lam * np.sech(100 * (-.5 + v[i]/w[i]))**2 * (3/4 + (1/4)*np.tanh(100 * (-.5 + v[i]/w[i])))]])

def F(self, i):
if self.env.params_dict["Flux"]=="UP":
return np.array([
self.w[i] + self.lam * self.thetas[i] * self.w[i] + self.lam * (self.env.f(self.u_up[i]) - self.thetas[i-1]*self.w[i-1] - self.env.f(self.u_up[i-1])),
self.v[i] + self.lam * .5 * self.thetas[i]**2 * self.w[i] + self.lam * .5 * (self.env.f(self.u_up[i]) - self.thetas[i-1]**2 * self.w[i-1] - self.env.f(self.u_up[i-1]))
])

elif self.env.params_dict["Flux"]=="LF":
return np.array([
self.w[i] + self.lam * (self.thetas[i+1] * (self.f(self.w[i+1] + self.u_up[i+1]) - self.f(i+1) - self.alpha*self.w[i+1])
+ 2*self.alpha*self.thetas[i]*self.w[i] - self.thetas[i-1] * (self.f(self.w[i-1] + self.u_up[i-1]) - self.f(i-1) + self.alpha*self.w[i-1]))
+ self.lam * ((self.f(i+1) - self.f(i-1)) - self.alpha*(self.u_up[i+1] - 2*self.u_up[i] + self.u_up[i-1])),

self.v[i] + self.lam * ((.5*self.thetas[i+1]**2) * (self.f(self.w[i+1] + self.u_up[i+1]) - self.f(i+1) - self.alpha*self.w[i+1])
+ 2*self.alpha*(.5*self.thetas[i]**2)*self.w[i] - (.5*self.thetas[i-1]**2) * (self.f(self.w[i-1] + self.u_up[i-1]) - self.f(i-1) + self.alpha*self.w[i-1]))
+ .5*self.lam * ((self.f(i+1) - self.f(i-1)) - self.alpha*(self.u_up[i+1] - 2*self.u_up[i] + self.u_up[i-1]))
])

def compute_J(self, i, epsilon=1e-6): #Compute the Jacobian
if self.env.params_dict["Theta_choice_method"] == "MinMax":
if self.env.params_dict["Flux"]=="UP":
if self.env.params_dict["f"] == "linear_advection":
if np.abs(self.w[i])>epsilon :
if self.v[i]/self.w[i] > self.theta_min:
J = np.array([[1, self.alpha * self.lam], [-(self.v[i]**2)*self.alpha*self.lam / (2*(self.w[i]**2)), 1 + self.alpha*self.lam/self.w[i]]])
else :
J = np.array([[1+self.lam*self.theta_min*self.alpha, 0], [self.lam * self.alpha*(self.theta_min**2) / 2, 1]])
else :
J = np.array([[1+self.lam*self.theta_st*self.alpha, 0], [self.lam * (self.theta_st**2) / 2, 1]])

elif self.env.params_dict["Flux"] == "LF":
pass


if self.env.params_dict["Theta_choice_method"] == "Smoothing":
x = self.tanh_factor
J = np.array([[1 - 25 * self.v[i] * self.lam * self.sech(100*(-.5+ self.v[i]/self.w[i]))**2 / self.w[i] + self.lam * (.75 + .25*np.tanh(x*(-.5 + self.v[i]/self.w[i]))) ,
25 * self.lam * self.sech(100*(-.5 + self.v[i]/self.w[i]))**2],
[-(25 * self.v[i] * self.lam * self.sech(100*(-.5 + self.v[i]/self.w[i]))**2 * (.75 + .25*np.tanh(x*(-.5 + self.v[i]/self.w[i])))) / self.w[i] + .5 * self.lam * (.75 + .25*np.tanh(x * (-.5 + self.v[i]/self.w[i])))**2 ,
1 + 25 * self.lam * self.sech(100 * (-.5 + self.v[i]/self.w[i]))**2 * (3/4 + (1/4)*np.tanh(x * (-.5 + self.v[i]/self.w[i])))]])
#print(np.linalg.norm(J))
return J


def Newton(self, epsilon=1e-6, maxiter=5e2):
def Newton(self, epsilon=1e-6, maxiter=10):
w_ = np.copy(self.u_down)
w_[0] = 0
#v = np.empty_like(self.u_down)
v = self.u_til - self.u_down
v[0] = 0 #
w = np.zeros_like(w_)
self.v = self.u_til - self.u_down
self.v[0] = 0 #
self.w = np.zeros_like(w_)
#(Neumann at the left)
iter=0

while np.linalg.norm(w - w_) > epsilon and maxiter>iter:
while np.linalg.norm(self.w - w_) > epsilon and maxiter>iter:
for i in range(1,self.env.mesh.Nx +1):
w_[i] = w[i]
w_[i] = self.w[i]

near0_eps = 1e-12
#if w[i] != 0:
if np.abs(self.w[i]) > near0_eps:
#update thetas
self.theta_choice(i, epsilon=self.env.params_dict["Theta_choice_epsilon"])
#compute step (linear system resolution)
s = np.linalg.solve(self.compute_J(i, epsilon=self.env.params_dict["Theta_choice_epsilon"]), -self.F(w,v,i))
else:
self.thetas[i] = 1.
s = np.linalg.solve(np.array([[1+self.lam,0],[self.lam/2,1]]), -self.F(i))

#update thetas
self.theta_choice(v[i],w[i],i, epsilon=self.env.default_dict["Theta_choice_epsilon"])
#compute step (linear system resolution)
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]
self.w[i] += s[0]
self.v[i] += s[1]
#print(np.mean(w))
iter += 1

if self.env.default_dict["print_Newton_iter"]==True:
if self.env.params_dict["print_Newton_iter"]==True:
print(iter, end=" ")

for i in range(w.size):
self.theta_choice(v[i],w[i],i, epsilon=self.env.default_dict["Theta_choice_epsilon"])
self.u_til = v

for i in range(self.w.size):
self.theta_choice(i, epsilon=self.env.params_dict["Theta_choice_epsilon"])

"""
def fixed_point_iter(self, theta_left, u_down_left, u_down,
w_left, i, epsilon=1e-6):
w = w_left #
w_ = epsilon + w #Just to initialize
theta = self.theta_st #Try wit hthe previous time step theta
while np.abs(w_ - w) >= epsilon:
w_ = w
w = (u_down - theta_left*w_left - (u_down_left + w_left)) / (1 + self.lam*theta)
v = (-self.lam/2) * (u_down - theta_left**2 * w_left - (u_down_left + w_left) + theta**2*w)
self.theta_choice(v,w,i, epsilon=self.env.default_dict["Theta_choice_epsilon"])
theta = self.thetas[i]
return theta, w
"""


def update_thetas(self):

if self.theta_computation == "Newton":
self.Newton(epsilon=self.env.default_dict["Method_convergence_epsilon"])

"""
elif self.theta_computation == "Fixed-Point":
#self.u_down[0] = self.left_bound_val
#w_left = self.u_up[0] - self.u_down[0] #Boundary condition: the value at the left is constant
w_left = 0 #For now, we work with a constant flux at the left, so the first w_left is equal to 0
self.Newton(epsilon=self.env.params_dict["Method_convergence_epsilon"])

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, i, epsilon=self.env.default_dict["Method_convergence_epsilon"])
"""
else:
print("Wrong theta choice method type")


def SATh_Scheme(self):

t = 0
self.thetas = np.ones(self.env.mesh.Nx+1) * self.theta_st
self.u_down = self.env.funcs.init_sol.copy()
self.u_up = np.empty_like(self.u_down)
self.u_up[0] = self.u_down[0]

if self.env.default_dict["Timer"]==True:
if self.env.params_dict["Timer"]==True:
self.timer("set")

while (t<self.env.tf):
Expand All @@ -228,10 +247,10 @@ def SATh_Scheme(self):
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)
self.u_til = self.v + self.u_up #update u_til

t += self.env.mesh.dt
if self.env.default_dict["Timer"]==True:
if self.env.params_dict["Timer"]==True:
self.timer()

return self.u_up
Expand Down
3 changes: 2 additions & 1 deletion args.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,9 +11,10 @@
"Theta_min":0,
"Theta_max":1,
"Scheme": "SATh",
"Theta_computation": "Fixed-Point",
"Theta_computation": "Newton",
"Theta_choice_epsilon":1e-6,
"Method_convergence_epsilon":1e-6,
"Theta_choice_method": "MinMax",
"Timer":0,
"print_Newton_iter":0,
"Path": "pictures",
Expand Down

0 comments on commit faf8db3

Please sign in to comment.