Add files via upload

eigenvalue_problem_ocean_basins.py

@@ -0,0 +1,94 @@
+import numpy as np 
+import matplotlib.pyplot as plt
+
+"""
+This is the code that solves the ocean basins eigenproblem. This is done with
+multiple shifts.
+"""
+
+plt.rcParams.update({'font.size': 30})
+
+from firedrake import *
+Lx, Ly = 1, 1
+n0 = 50
+mesh = RectangleMesh(n0, n0, Lx, Ly, reorder=None)
+V = FunctionSpace(mesh, "CG", 1)
+bc = DirichletBC(V, 0, "on_boundary")
+n = 300
+phi, psi = TestFunction(V), TrialFunction(V)
+
+opts = {"eps_gen_non_hermitian": None, "eps_largest_imaginary": None,
+        "st_type": "shift", "eps_target": None,
+        "st_pc_factor_shift_type": "NONZERO"} # options for solver
+
+eigenproblem_2_shift = LinearEigenproblem(
+        A=phi*psi.dx(0)*dx,
+        M=-inner(grad(psi), grad(phi))*dx - psi*phi*dx,
+        bcs=bc, bc_shift=2, restrict=False)
+eigensolver_2_shift = LinearEigensolver(eigenproblem_2_shift, n_evals=n,
+                                        solver_parameters=opts) 
+nconv_2_shift = eigensolver_2_shift.solve()
+
+
+eigenproblem_0_shift = LinearEigenproblem(
+        A=phi*psi.dx(0)*dx,
+        M=-inner(grad(psi), grad(phi))*dx - psi*phi*dx,
+        bcs=bc, bc_shift=0.0, restrict=False)
+
+eigensolver_0_shift = LinearEigensolver(eigenproblem_0_shift, n_evals=n,
+                                        solver_parameters=opts) 
+nconv_0_shift = eigensolver_0_shift.solve()
+
+print("Number of eigenvalues converged (shift=2):", nconv_2_shift)
+print("Number of eigenvalues converged (shift=0):", nconv_0_shift)
+print("Number of eigenvalues searched for:", n)
+
+min_n_conv = min(nconv_2_shift, nconv_0_shift)
+evals_2_shift_real, evals_2_shift_imag = np.zeros(min_n_conv, dtype=complex), np.zeros(min_n_conv, dtype=complex)
+evals_0_shift_real, evals_0_shift_imag = np.zeros(min_n_conv, dtype=complex), np.zeros(min_n_conv, dtype=complex)
+for i in range(min_n_conv):
+    eval = eigensolver_2_shift.eigenvalue(i)
+    evals_2_shift_real[i] = np.real(eval)
+    evals_2_shift_imag[i] = np.imag(eval)
+    eval2 = eigensolver_0_shift.eigenvalue(i)
+    evals_0_shift_real[i] = np.real(eval2)
+    evals_0_shift_imag[i] = np.imag(eval2)
+
+plt.scatter(evals_2_shift_real, evals_2_shift_imag, color="b", label=r"Shift=2")
+plt.scatter(evals_0_shift_real, evals_0_shift_imag, color="r", label=r"Shift=0")
+plt.xscale("log")
+plt.title(rf"First {min_n_conv} Eigenvalues [Sorted By Largest Imaginary Part], Created Using Different Shift Values", wrap=True)
+plt.xlabel("Real Part")
+plt.ylabel("Imaginary Part")
+plt.legend()
+plt.show()
+
+
+eigenproblem_res = LinearEigenproblem(
+        A=phi*psi.dx(0)*dx,
+        M=-inner(grad(psi), grad(phi))*dx - psi*phi*dx,
+        bcs=bc, bc_shift=2, restrict=True)
+eigensolver_res = LinearEigensolver(eigenproblem_res, n_evals=n, solver_parameters=opts)
+nconv_res = eigensolver_res.solve()
+
+print("Number of eigenvalues converged (restricted):", nconv_res)
+
+evals_real_res = np.zeros(nconv_res, dtype=complex)
+evals_imag_res = np.zeros(nconv_res, dtype=complex)
+for i in range(nconv_res):
+    eval = eigensolver_res.eigenvalue(i)
+    evals_real_res[i] = np.real(eval)
+    evals_imag_res[i] = np.imag(eval)
+
+overall_min_nconv = min(min_n_conv, nconv_res)
+
+# plotting just up to the min of all three converged values?
+plt.scatter(evals_2_shift_real[:overall_min_nconv], evals_2_shift_imag[:overall_min_nconv], color="b", label=r"Shift=2")
+plt.scatter(evals_0_shift_real[:overall_min_nconv], evals_0_shift_imag[:overall_min_nconv], color="r", label=r"Shift=0")
+plt.scatter(evals_real_res[:overall_min_nconv], evals_imag_res[:overall_min_nconv], color="g", label=r"Restricted")
+plt.xscale("log")
+plt.title(rf"First {overall_min_nconv} Eigenvalues for the Problem, in the Restricted and Unrestricted Eigensolver", wrap=True)
+plt.xlabel("Real Part")
+plt.ylabel("Imaginary Part")
+plt.legend()
+plt.show()

manufactured_solution.py

@@ -0,0 +1,66 @@
+"""
+This is to solve the problem, with a manufactured solution, in the comparison
+to function space section.
+"""
+from firedrake import *
+import matplotlib.pyplot as plt
+from firedrake.pyplot import tripcolor, tricontour
+import pygments
+plt.rcParams.update({'font.size': 26})
+
+mesh = UnitSquareMesh(16, 16)
+x, y = SpatialCoordinate(mesh)
+
+# solving in the unrestricted case
+V = FunctionSpace(mesh, "CG", 2)
+u, v = TrialFunction(V), TestFunction(V)
+f = Function(V).interpolate(-2 * (y**3 - 1.5 * y**2) + (x-x**2) * (6*y -3))
+G = inner(f, v) * dx
+a = inner(grad(u), grad(v)) * dx
+bc_left = DirichletBC(V, 0.0, 1)
+bc_right = DirichletBC(V, 0.0, 2)
+
+u = Function(V)
+solve(a == G, u, bcs=[bc_left, bc_right], restrict=False) 
+
+V_res = RestrictedFunctionSpace(V, boundary_set=[1, 2])
+u_res, v_res = TrialFunction(V_res), TestFunction(V_res)
+f_res = Function(V_res).interpolate(-2 * (y**3 - 1.5 * y**2) + (x-x**2) * (6*y -3))
+G_res = inner(f_res, v_res) * dx
+a_res = inner(grad(u_res), grad(v_res)) * dx
+
+bc_left_res = DirichletBC(V_res, 0.0, 1)
+bc_right_res = DirichletBC(V_res, 0.0, 2)
+
+u_res = Function(V_res)
+solve(a_res == G_res, u_res, bcs=[bc_left_res, bc_right_res], restrict=False)
+
+# interpolating the known solution into the function space. 
+u_solution = Function(V).interpolate(-1 * x * (1-x) * (y**3 - 1.5 * y**2))
+
+# printing norms of interest. 
+print("Error norm between restricted and unrestricted solutions:", errornorm(u, u_res)) 
+print("Error norm between unrestricted and true solutions:", errornorm(u, u_solution))
+print("Error norm between restricted and true solutions:", errornorm(u_res, u_solution))
+
+# plotting
+fig, [ax1, ax2, ax3] = plt.subplots(nrows=1, ncols=3)
+plot_u = tripcolor(u, axes=ax1)
+plot_solution = tripcolor(u_solution, axes=ax2)
+plot_ures = tripcolor(u_res, axes=ax3)
+
+ax1.set_xlabel("x")
+ax1.set_ylabel("y")
+ax2.set_xlabel("x")
+ax2.set_ylabel("y")
+ax2.yaxis.set_label_coords(-0.05, 0.5)
+ax3.set_xlabel("x")
+ax3.set_ylabel("y")
+ax3.yaxis.set_label_coords(-0.05, 0.5)
+ax1.set_title("Numerical Solution \n(Unrestricted)", wrap=True)
+ax2.set_title("Interpolation of True \nSolution", wrap=True)
+ax3.set_title("Numerical Solution \n(Restricted)", wrap=True)
+fig.subplots_adjust(right=0.8)
+cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
+fig.colorbar(plot_ures, cax=cbar_ax)
+plt.show()

poisson_eigenvalue_problem.py

@@ -0,0 +1,67 @@
+"""
+This is basically a clone of the eigensolver test in Firedrake, with images
+added.
+"""
+from firedrake import * 
+import matplotlib.pyplot as plt
+import numpy as np
+plt.rcParams.update({'font.size': 22})
+
+n = 10
+degree = 4
+
+mesh = IntervalMesh(n, 0, pi) # setting the mesh between 0 and pi
+V = FunctionSpace(mesh, "CG", degree)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+a = (inner(grad(u), grad(v))) * dx
+
+bc = DirichletBC(V, 0.0, "on_boundary") # this is equivalent to a bc on subdomain 1 and 2
+
+# the default M is the one in this problem, so no need to specify in construction
+eigenprob_res = LinearEigenproblem(a, bcs=bc, bc_shift=1, restrict=True)
+eigenprob = LinearEigenproblem(a, bcs=bc, bc_shift=70, restrict=False)
+
+n_evals = 10
+
+eigensolver_res = LinearEigensolver(
+    eigenprob_res, n_evals, solver_parameters={"eps_largest_real": None}
+)
+nconv_res = eigensolver_res.solve()
+
+eigensolver = LinearEigensolver(
+    eigenprob, n_evals, solver_parameters={"eps_largest_real": None}
+)
+nconv = eigensolver.solve()
+
+print("Number of eigenvalues found (restricted):", nconv_res)
+print("Number of eigenvalues found (unrestricted):", nconv)
+
+estimates_res = np.zeros(nconv_res, dtype=complex)
+estimates = np.zeros(nconv, dtype=complex)
+for k in range(max(nconv_res, nconv)):
+    if k < nconv_res:
+        estimates_res[k] = eigensolver_res.eigenvalue(k)
+    if k < nconv:
+        estimates[k] = eigensolver.eigenvalue(k)
+
+eigenmode_real, eigenmode_imag = eigensolver.eigenfunction(0)
+eigenmode_r_res, eigenmode_i_res = eigensolver_res.eigenfunction(0)
+
+true_values = [x**2 for x in range(1, nconv_res+1)] # known eigenvalues
+
+for eval in true_values[:n_evals]:
+    plt.vlines(eval, ymin=-0.001, ymax=0.001, linestyle="--", color="black", zorder=-1)
+plt.scatter(np.real(estimates_res)[:n_evals], 0.001 * np.ones(nconv_res)[:n_evals], 
+            s=200, color="b", marker="x", label="Restricted", zorder=1)
+plt.scatter(np.real(estimates)[:n_evals], -0.001 * np.ones(nconv)[:n_evals], 
+            s=200, color="r", marker="s", label="Unrestricted", zorder=1)
+plt.scatter(np.real(true_values)[:n_evals], np.zeros(len(true_values))[:n_evals], 
+            s=200, color="g", marker="*", label="Exact", zorder=1)
+
+plt.title(f"First {n_evals} Converged Eigenvalues for the 1D Poisson Eigenproblem")
+plt.tick_params(axis='y',  which='both', left=False,  right=False, labelleft=False) 
+plt.xlabel(r"$\lambda$ (real)")
+plt.legend(loc=(0.79, 0.2))
+plt.show()

solve_small_poisson_problem.py

@@ -0,0 +1,28 @@
+"""
+This code is to demonstrate the use of Firedrake in solving a small Poisson 
+equation
+"""
+
+from firedrake import *
+mesh = UnitSquareMesh(1, 1)
+x, y = SpatialCoordinate(mesh)
+# create Lagrange element of degree 4
+V = FunctionSpace(mesh, "Lagrange", 4)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+
+# 1 indicates the left side of mesh, 2 for right side
+bc_left = DirichletBC(V, 0.0, 1)
+bc_right = DirichletBC(V, 0.0, 2)
+
+# create bilinear form, dx means integral
+a = inner(grad(u), grad(v)) * dx
+
+# create linear form
+f = Function(V).interpolate(x**2)
+G = inner(f, v) * dx 
+
+# solve, u_sol holds solution
+u_sol = Function(V)
+solve(a == G, u_sol, bcs=[bc_left, bc_right])

solve_small_poisson_problem_pres.py

@@ -0,0 +1,22 @@
+"""
+Code with comments/blank lines removed
+"""
+
+from firedrake import *
+mesh = UnitSquareMesh(1, 1)
+x, y = SpatialCoordinate(mesh)
+V = FunctionSpace(mesh, "Lagrange", 2)
+u = TrialFunction(V)
+v = TestFunction(V)
+bc_left = DirichletBC(V, 0.0, 1)
+bc_right = DirichletBC(V, 0.0, 2)
+bcs=[bc_left, bc_right]
+a = inner(grad(u), grad(v)) * dx
+f = Function(V).interpolate(x**2)
+G = inner(f, v) * dx 
+K = assemble(a, bcs=bcs)
+u_sol = Function(V)
+solve(a == G, u_sol, bcs=bcs)
+
+print(K.M.values)
+print("u values:", u_sol.dat.data)