#==================
# solver.py (with adaptive alpha strategy)
#==================
import numpy as np
from scipy.sparse.linalg import spsolve
import time
class NonlinearSolver:
"""Incremental-iterative solver with adaptive alpha strategy"""
def __init__(self, assembly, n_increments=10, tol=1e-6, max_iter=20,
ls_max=10, ls_alpha=0.5, ls_c1=1e-4,
alpha_reduction=0.5, max_alpha_reductions=3):
self.assembly = assembly
self.n_increments = n_increments
self.tol = tol
self.max_iter = max_iter
self.ls_max = ls_max
self.ls_alpha = ls_alpha
self.ls_c1 = ls_c1
self.alpha_reduction = alpha_reduction
self.max_alpha_reductions = max_alpha_reductions
self.u = np.zeros(assembly.n_dof)
self.residual_norms = []
self.convergence_data = []
def solve(self):
"""Solve with adaptive alpha strategy"""
# Initialize displacement vector
self.assembly.u = np.zeros(self.assembly.n_dof)
# Apply gravity (constant)
self.assembly.F_ext = np.zeros(self.assembly.n_dof)
self.assembly.apply_gravity()
# Initialize load factor
load_factor = 0.0
load_increment = 1.0 / self.n_increments
# Time the solution
start_time = time.time()
for inc in range(1, self.n_increments + 1):
# Update load factor
load_factor += load_increment
print(f"\n--- Load Increment {inc}/{self.n_increments}
(λ={load_factor:.2f}) ---")
# Apply fluid pressures
self.assembly.apply_hydrostatic_pressure(load_factor * 95)
self.assembly.apply_uplift_pressure(
upstream_level=load_factor * 1000*9.81*95,
downstream_level=0
)
# Save state at start of increment
current_state = self.assembly.save_state()
alpha_reductions = 0
converged = False
while alpha_reductions <= self.max_alpha_reductions and not converged:
# Restore to start of increment
self.assembly.restore_state(current_state)
# Adjust line search parameters
current_ls_alpha = self.ls_alpha * (self.alpha_reduction **
alpha_reductions)
# Attempt to solve increment with current alpha
converged = self.solve_increment(current_ls_alpha)
if not converged:
alpha_reductions += 1
print(f"Alpha reduction #{alpha_reductions} (new
α={current_ls_alpha:.4f})")
if converged:
print(f"Converged with α={current_ls_alpha:.4f}")
self.convergence_data.append((load_factor, alpha_reductions))
else:
print(f"Warning: Failed to converge in increment {inc}")
break
total_time = time.time() - start_time
print(f"\nSolution completed in {total_time:.2f} seconds")
print(f"Total iterations: {len(self.residual_norms)}")
print(f"Final load factor: {load_factor:.4f}")
return self.u
def solve_increment(self, current_ls_alpha):
"""Solve current increment with specified line search alpha"""
# Initial stiffness assembly
self.assembly.assemble_stiffness(tangent=False)
for iter in range(self.max_iter):
# Compute internal forces
F_int = self.assembly.compute_internal_forces()
# Compute residual
R = self.assembly.F_ext - F_int
norm_R = np.linalg.norm(R)
self.residual_norms.append(norm_R)
# Apply boundary conditions
K_mod, R_mod = self.assembly.apply_boundary_conditions(self.assembly.K,
R)
# Solve for displacement increment
try:
du = spsolve(K_mod, R_mod)
except Exception as e:
print(f"Linear solver error: {str(e)}")
return False
# Perform line search with current alpha
alpha = self.armijo_line_search(du, R, current_ls_alpha)
# Update displacements
self.assembly.u += alpha * du
self.u = self.assembly.u
# Check convergence
norm_du = np.linalg.norm(alpha * du)
print(f"Iter {iter+1}: |R| = {norm_R:.4e}, |Δu| = {norm_du:.4e}, α =
{alpha:.4f}")
if norm_R < self.tol and norm_du < self.tol:
return True
# Update tangent stiffness
self.assembly.assemble_stiffness(tangent=True)
return False
def armijo_line_search(self, du, R0, base_alpha):
"""Armijo line search with specified base alpha"""
alpha = 1.0 # Start with full step
u_prev = self.assembly.u.copy()
norm_R0 = np.linalg.norm(R0)
min_alpha = 0.01 # Minimum allowed step size
for i in range(self.ls_max):
# Trial displacement
self.assembly.u = u_prev + alpha * du
# Compute new residual
F_int = self.assembly.compute_internal_forces()
R_new = self.assembly.F_ext - F_int
norm_R_new = np.linalg.norm(R_new)
# Armijo condition: sufficient decrease
if norm_R_new <= norm_R0 + self.ls_c1 * alpha * np.dot(R0, du):
return alpha
# Reduce step size using current base reduction factor
alpha *= base_alpha
if alpha < min_alpha:
return min_alpha
return alpha
